“Calhoun 


Institutional Archive of the Naval Postgraduate School 





Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


2001-06 


Angular rate estimation for multi-body 
Spacecraft attitude control 


Palermo, William J. 


Monterey, California. Naval Postgraduate School 
http://hdl.handle.net/10945/2389 


This publication is a work of the U.S. Government as defined in Title 17, United 
States Code, Section 101. Copyright protection is not available for this work in the 
United States. 


Downloaded from NPS Archive: Calhoun 


Calhoun is the Naval Postgraduate School's public access digital repository for 


\§ D U DL EY research materials and institutional publications created by the NPS community. 
«iis Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


NNN KNOX appointed -- and published -- scholarly author. 


LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 


http://www.nps.edu/library Monterey, California USA 93943 


NAVAL POSTGRADUATE SCHOOL 
Monterey, California 





THESIS 


ANGULAR RATE ESTIMATION FOR MULTI-BODY 
SPACECRAFT ATTITUDE CONTROL 


by 
William J. Palermo 
June 2001 


Thesis Advisor: Brij N. Agrawal 
Second Reader: Harold A Titus 








Approved for public release; distribution is unlimited 


REPORT DOCUMENTATION PAGE 


Public reporting burden for this collection of information is estimated to average | hour per response, including 
the time for reviewing instruction, searching existing data sources, gathering and maintaining the data needed, and 
completing and reviewing the collection of information. Send comments regarding this burden estimate or any 
other aspect of this collection of information, including suggestions for reducing this burden, to Washington 
headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 
1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project 
(0704-0188) Washington DC 20503. 


1. AGENCY USE ONLY (Leave blank) | 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 
June 2001 Engineer’s Thesis 
4. TITLE AND SUBTITLE: Angular Rate Estimation for Multi-Body Spacecraft | 5. FUNDING NUMBERS 
Attitude Control 


6. AUTHOR(S) William J. Palermo 


7, PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING 
Naval Postgraduate School ORGANIZATION REPORT 
Monterey, CA 93943-5000 NUMBER 

9. SPONSORING / MONITORING AGENCY NAME(S) AND ADDRESS(ES) | 10. SPONSORING / MONITORING 
N/A AGENCY REPORT NUMBER 


11. SUPPLEMENTARY NOTES The views expressed in this thesis are those of the author and do not reflect the official 
policy or position of the Department of Defense or the U.S. Government. 


12a. DISTRIBUTION / AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE 
Approved for public release; distribution is unlimited 


13. ABSTRACT (maximum 200 words) 


Spacecraft with high performance attitude control systems requirements have traditionally relied on imperfect mechanical 
gyroscopes for primary attitude determination. Gyro bias errors are corrected with a Kalman filter algorithm that uses updates 
from precise attitude sensors like star trackers. Gyroscopes, however, have a tendency to degrade or fail on orbit, becoming a 
life limiting factor for many satellites. When errors become erratic, pointing accuracy may be lost during short star gaps. 
Unpredictable gyro degradations have impacted NASA spacecraft missions such as Skylab and Hubble Space Telescope as 
well as several DoD and ESA satellites. An alternative source of angular rate information is a software implemented real time 
dynamic model. Inputs to the model from internal sensors and known spacecraft parameters enable the tracking of total system 
angular momentum from which body rates can be determined. With this technique, the Kalman filter algorithm provides error 
corrections to the dynamic model. The accuracy of internal sensors and input parameters determine the effectiveness of this 
angular rate estimation technique. This thesis presents the background for understanding and implementation of this technique 
into a representative attitude determination system. The system is incorporated into an attitude simulation model cveloped in 
SIMULINK to evaluate the effects of dynamic modeling errors and sensor inaccuracies. Results are presented that indicate that 
real time dynamic modeling is an effective method of angular rate determination for maneuvering multi-body spacecraft 
attitude control systems. 


14. SUBJECT TERMS dynamic gyro, Kalman filter, attitude determination, rate estimation, star | 15. NUMBER OF 
trackers, attitude simulation, multi-body dynamics, quaternions, MATLAB, SIMULINK PAGES 


16. PRICE CODE 


17. SECURITY 18. SECURITY 19. SECURITY 20. LIMITATION 
CLASSIFICATION OF CLASSIFICATION OF THIS CLASSIFICATION OF OF ABSTRACT 
REPORT PAGE ABSTRACT 
Unclassified Unclassified Unclassified UL 
NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89) 
Prescribed by ANSI Std. 239-18 





THIS PAGE INTENTIONALLY LEFT BLANK 


li 


Approved for public release; distribution is unlimited 


ANGULAR RATE ESTIMATION FOR MULTI-BODY SPACECRAFT 
ATTITUDE CONTROL 


William J. Palermo 
Lieutenant, United States Navy 
B.S., United States Naval Academy, 1992 


Submitted in partial fulfillment of the 
requirements for the degree of 


MASTER OF SCIENCE IN ASTRONAUTICAL ENGINEERING 
and 
AERONAUTICAL AND ASTRONAUTICAL ENGINEER 


from the 


NAVAL POSTGRADUATE SCHOOL 
June 2001 


Author: 





William J. Palermo 


Approved by: 





Brij N. Agrawal, Thesis Advisor 





Harold A. Titus, Second Reader 





Max F. Platzer, Chairman 
Department of Aeronautics and Astronautics 


iil 


THIS PAGE INTENTIONALLY LEFT BLANK 


iv 


ABSTRACT 


Spacecraft with high performance attitude control systems requirements have 
traditionally relied on imperfect mechanical gyroscopes for primary attitude 
determination. Gyro bias errors are corrected with a Kalman filter algorithm that uses 
updates from precise attitude sensors like star trackers. Gyroscopes, however, have a 
tendency to degrade or fail on orbit, becoming a life limiting factor for many satellites. 
When errors become erratic, pointing accuracy may be lost during short star gaps. 
Unpredictable gyro degradations have impacted NASA spacecraft missions such as 
Skylab and Hubble Space Telescope as well as several DoD and ESA satellites. An 
alternative source of angular rate information is a software implemented real time 
dynamic model. Inputs to the model from internal sensors and known spacecraft 
parameters enable the tracking of total system angular momentum from which body rates 
can be determined. With this technique, the Kalman filter algorithm provides error 
corrections to the dynamic model. The accuracy of internal sensors and input parameters 
determine the effectiveness of this angular rate estimation technique. This thesis presents 
the background for understanding and implementation of this technique into a 
representative attitude determination system. The system is incorporated into an attitude 
simulation model developed in SIMULINK to evaluate the effects of dynamic modeling 
errors and sensor inaccuracies. Results are presented that indicate that real time dynamic 
modeling is an effective method of angular rate determination for maneuvering multi- 


body spacecraft attitude control systems. 
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I. INTRODUCTION 


A. ANGULAR RATE ESTIMATION 


The attitude control systems of today’s high tech satellites require accurate 
angular rate knowledge for attitude propagation and control loop damping in order to 
meet pointing and tracking requirements. Gyroscopes offer varying degrees of precision 
for direct measurement of angular rates and have been the primary attitude determination 
sensor used by spacecraft for many years. However, high cost and low reliability has 
lead users and designers to explore other options. Using dynamic calculations, 
uncertainty based estimation algorithms or a combination of these two methods, onboard 
processors can estimate spacecraft angular rates without measuring them directly. 
Angular rate estimation techniques can be a viable alternative for back-up or even 


primary attitude determination system in many control schemes. 


B. OPTIMAL ESTIMATION 


“An optimal estimator is a computational algorithm that processes measurements 
to deduce a minimum error estimate, in accordance with stated optimality criterion, of the 
state of a system by utilizing: knowledge of system measurement dynamics, assumed 
Statistics of system noises and measurement errors, and initial condition information.” 
[Ref. 1] Estimation techniques provide filtering, smoothing and prediction of state 
variables in an imperfect model based on imperfect measurement update data. The most 
common optimal estimator used in stochastic systems is the Kalman filter. The dynamic 
model error and measurement error are assumed as zero mean Guassian white noise 
processes with known covariance. Estimators are commonly used even in systems where 


all state variables required by controller can be measured. 
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Figure 1. Basic Kalman Filter Block Diagram 


OBSERVATION 








Maybeck identifies three basic reasons why deterministic methods are insufficient 
for describing real systems. First, no mathematical model is perfect. There are always 
effects that are necessarily neglected to make the model practical and even modeled 
effects are only approximations to what is physically occurring. | Second, dynamic 
systems are driven not only by control inputs but also by disturbances which we can 
neither control nor model deterministically. Lastly sensors do not provide perfect and 


complete data. [Ref. 2] 


An important advantage of using an optimal estimator such as the Kalman filter is 
that the attitude determination output does not affect the controller design. Therefore the 
development of an optimal controller can be accomplished independently. An attitude 
control algorithm optimized for ideal state inputs will remain optimized with state 


estimates provided by an optimal estimator. 


c. ANGULAR RATE ESTIMATION FROM VECTOR MEASUREMENTS 


The concept of using attitude sensor data to produce an estimate of spacecraft 
angular rate without gyroscopes is not new. Gyroless attitude and angular rate estimation 
has been a prime concern for small inexpensive spacecraft that do not carry gyroscopes 
but still require rate information for attitude propagation and control. Estimation 
techniques also provide options for complex spacecraft that require back-up control 


modes in the event of gyro failures. The problem of angular rate estimation can be 


treated separately from attitude determination. Several reliable algorithms have been 


developed that produce angular rate using on-board processors. 


It is possible to extract angular rate directly from time derivatives of measured 
vectors resolved in the body coordinate frame and known in an inertial frame from a 
model or almanac data. This technique, however, requires at least two non-parallel 
vector measurements. It also exhibits time lag and produces very noisy data since the 
algorithms depend on derivatives of already noisy measurements. A better estimate of 
angular rate can be achieved by treating the problem as stochastic. For this method an 
estimating filter is applied that uses dynamic equations of motion to take advantage of 
past data in a recursive sense. This method has the additional advantage of being able to 
update rate estimates at time steps when only one vector measurement is available. The 


vector derivative is treated as a noisy measurement update to a Kalman rate estimator. 


A Kalman filter requires a linear state-space dynamic model but the dynamics of a 
spacecraft in the rotating body coordinate frame are coupled and nonlinear. Attitude 
determination schemes deal with this problem in different ways. Several methods have 
been proposed for interlacing two or three linear Kalman filters together to capture the 
nonlinear dynamics [Ref. 3]. It is shown that the coupled equations of motion can be 
written as a linear combination of the spacecraft angular rate and two other newly defined 
vectors whose components are nonlinear combinations of the angular rate vector 
elements. The differential equations for each of these new vectors can be written as 
linear combinations of the other and the angular rate. Adding white noise vectors to these 
new equations of motion produces a set of three stochastic linear models that can be used 
in separate Kalman filters. The filters are interlaced with their estimated outputs treated 


as deterministic inputs to each other. 


An effective method of dealing with non-linear dynamics is presented for use by 
the Pseudo-linear Kalman filter (PSELIKA) and the state-dependent algebraic Riccati 
equation filter (SDARE) [Ref. 4]. For these unique rate estimation techniques the 
equation of motion is transformed to express the nonlinear terms in angular rate as a 
product of a matrix whose elements depend on the components of the angular rate vector 


and the angular rate vector itself. The pseudolinear Kalman filter or PSELIKA is 
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implemented as an ordinary discrete Kalman filter with a time dependent state transition 
matrix whose angular rate dependent elements are formed from the current estimate of 
angular rate. The SDARE or state-dependent algebraic Riccati equation filter is also 
implemented as a discrete Kalman filter derived from the same representation of the 
dynamic equation. The Kalman gain used in this filter is computed from the solution of 
an algebraic Riccati equation involving the angular rate dependent matrix. This 
eliminates the need © propagate and update the filter state covariance which is normally 


used in the gain calculation. 


The pseudolinear filter concept can be extended to use quaternion measurements. 
In the attitude determination system on board the Rossi X-Ray Timing Explorer (RXTE) 
satellite advanced star tracker software directly provides information in the form of 
quaternions [Ref. 5]. The same propagation algorithm is used with measurement updates 
in the form of quaternions. Since the device yielding the attitude measurement also 
conducts a star search and identification process a time delay is introduced. Two 


algorithms are presented for overcoming this delay problem. 


D. FULL ORDER ESTIMATORS 


It is not necessary to treat the problem of angular rate estimation separately from 
attitude determination. In fact, attitude sensor data also requires filtering to smooth out 
measurement noise and produce a clean attitude estimates. High precision control laws 
generally require both attitude and angular rate information. These components make up 


the full system state vector. 


Full order estimators produce optimal state estimates in the case where imperfect 
measurement data related to all of the modeled states is continuously available. 
Spacecraft with a triad of rate gyroscopes and attitude sensors have the required state 
measurements but the bandwidth of the attitude sensors is normally too low for the 
control system. Additionally, there may be periodic gaps in sensor coverage and certain 
sensors do not provide reliable attitude information in all three orthogonal dimensions. 


Therefore standard full order state estimators are generally not used. 


In most practical implementations of gyro-based attitude control systems rate 
information is used in the propagation stage of an attitude estimator based on kinematic 
equations. The estimator utilizes discrete vector observations, resolved in both body- 
fixed coordinates and a reference frame, to estimate spacecraft attitude and gyro drift 


rates. 


In an ideal attitude determination system where perfect knowledge of the 
spacecraft angular rates is available, the accuracy of the kinematic model for determining 
attitude is limited only by the processor time step and the initial attitude state. Attitude 
sensors can be used to periodically reinitialize the kinematic model but the data they 
provide is discontinuous and corrupted by environmental effects and sensor design 


imperfections. 


In real attitude determination systems, the source of angular rate information may 
be noisy as well as biased, as is the case with spacecraft gyroscopes. This results in 
additional errors in the spacecraft attitude determination system. For control systems 
where attitude and angular rate information is critical, optimal estimation techniques are 
employed to combat the effects of model and sensor inaccuracies. As gyros degrade the 
estimator can be made more robust to plant error at the expense of attitude determination 
accuracy. At some point the gyro outputs may become too erratic to meet control system 


requirements. 


E. REDUCED ORDER ESTIMATORS 


Obtaining satellite angular rate estimates without the use of rate gyroscopes or 
other deterministic rate data can be accomplished with a reduced order estimator. This 
class of filters produces estimates of all modeled state variables when only a subset is 
directly related to the measurement data. In this case, attitude sensor information is 
available but no direct rate measurement is performed. Since gyro data is not available 


no update to the spacecraft rate is available between attitude measurements. 


The reduced estimator Kalman filter is formulated from the state space dynamic 
equations of motion linearized about the current estimate of state. Again, for spacecraft 


attitude control the state vector includes both the attitude and angular rate. The 


deterministic inputs to the model are the known externally applied moments usually 
consisting only of control torques. Unmodeled disturbances and other dynamic effects 
must be accounted for by a robust plant noise model. The estimates of attitude and 
angular rate are propagated forward by the dynamic model with discrete corrections from 


the attitude sensors. 


The performance of the reduced Kalman estimator can be improved by including 
higher order dynamic effects in the system plant model and including disturbance torques 
as deterministic inputs. However this costs extra processing power, requires additional 
sensors and cannot account for all of the unknowns. Additional dynamic complexities 
are introduced in multi-body satellites that have time varying moments of inertia, 
changing centers of mass and flexibility modes. Determining the magnitude and effects 
of disturbances and modeling simplifications for a particular spacecraft is an important 
design consideration that often requires rough calculations, simulation and engineering 
judgment. In general, since disturbances are of low amplitude and low frequency 


compared to control torques, their effects can be accounted for as plant error by the filter. 


One approach that avoids the use of the typically uncertain dynamic model 
altogether is to treat the spacecraft as a noncooperative target whose rotation must be 
tracked by a stochastic estimator [Ref. 6]. This concept is borrowed from tracking 
theory where it has been widely used to estimate target motion. This algorithm adds 
angular acceleration to the state vector and substitutes attitude states with the integrated 
rate parameters to formulate a nine state linear Kalman filter. Since this set of variables 
serves aS an approximate third-order attitude parameterization, the size must be 
controlled by sampling interval. Instead of using dynamics, time propagation of the 
estimated state variables is performed in the proposed filter by modeling the spacecraft 
angular acceleration as an exponentially autocorrelated stochastic process and using a 


polynomial kinematic model. 


In reduced order estimators, time steps for state propagation must be kept small to 
minimize the effects of dynamic simplifications and linearization. The primary limitation 
of the reduced estimator becomes evident during prolonged intervals between attitude 


sensor measurement updates. Gaps in attitude sensor coverage must be kept short enough 
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to ensure that estimator errors do not grow beyond required precision of the attitude 


control system. 


F. STAR SENSOR BASED ESTIMATORS 


Bandwidth and accuracy limitations of attitude sensors have precluded the use of 
reduced order estimators for spacecraft missions with high performance pointing and 
tracking requirements. Aside from star trackers, most sensors do not provide accurate 
and stable enough attitude data for high performance tracking requirements. 
Additionally, these sensors provide discrete measurements unlike gyroscopes which 
produce nearly continuous data and can be sampled at extremely high frequency limited 
only by processing capabilities. However, recent advancements in sensor technology 
suggest that star trackers can provide updates accurately and frequently enough to be 
considered for use as primary attitude sensors for a wide range of spacecraft missions. 
Key technological improvements include a wide field of view, high sensitivity, low noise 


equivalent angles and high bandwidth iteration rate. 


It has been proposed that advanced high-bandwidth star sensors can be used as the 
primary sensors for on-orbit attitude rate determination in place of angular rate 
gyroscopes [Ref. 7]. A Kalman filter algorithm, using updates from star tracker 
measurements, is presented that tracks errors in the attitude and angular rate for a 
kinematic based state space model. The spacecraft angular rate used for attitude 
propagation is separated into a nominal component and a component due to control and 
disturbance moments. The algorithm assumes that the bandwidths of the control system 
and disturbance effects are at least an order of magnitude smaller than the measurement 
bandwidth. This allows the rate of change in spacecraft angular rate due to unmodeled 
dynamic effects to be modeled as a first-order Gauss-Markov process. Nonlinear 
kinematics are used in the attitude propagation phase which is performed external to the 
estimator. The errors in the this model are then tracked using an ordinary linear discrete 


Kalman filter. 


G. ANALYTICAL RATE DETERMINATION 


Another method of obtaining estimates of spacecraft angular rates is through 
direct calculation from the dynamic equations of motion for the system. This rate 
calculation is performed in real time based on known, derived and sensed internal 
parameters of the spacecraft. The software implemented dynamic model can be adjusted 
to include varying levels of complexity for multi-body spacecraft that operate in different 
configurations. The Aerospace Corporation has patented one methodology for this 


process called the “Pseudo Gyro” [Ref. 8]. 


Information from internal sensors that detect relative orientations and angular 
rates of momentum exchange devices and appendages are used to continually update the 
parameters used in the dynamic calculations. These parameters determine component 
contributions to the system angular momentum and inertia dyadic. Control torques and 
modeled disturbances are integrated to capture external dynamic effects. Using the total 
system angular momentum and inertia dyadic and the relative momentum of internal 
components and appendages from the dynamic model the spacecraft angular rate can be 
calculated. If the appendage masses are small or their relative motions are slow the 
spacecraft inertia matrix is slow to change. The moment of inertia calculations can be 
performed at a lower frequency to save processing time. Cross product effects due to the 
rotating coordinate frame must also be accounted for in the dynamic model since all 


measurements are taken in the body frame. 


The accuracy of the calculated rate is dependent on the quality of the dynamic 
model and the sensor information available. Error sources include imperfect knowledge 
of component or appendage inertia matricies and mass centers as well as relative angular 
position and rate data from the internal sensors. Often relative rates are not measured 
directly but derived from position encoders which adds extra noise to the momentum 
along the axis of rotation of the appendage. Errors are also introduced through external 
disturbance effects since they cannot be perfectly modeled. All of these dynamic 
influences effect the model precision in varying degrees. Additionally, due to the 
finite/discrete processing capabilities of the computer performing the calculations there 


will be a slight drift over time from the true state even in the case of perfect sensors and 


input data. This drift is obviously minimized by increasing the frequency of the discrete 


model calculations. 


H. ATTITUDE ESTIMATION FROM CALCULATED RATE 


The rate estimate produced by the dynamic calculations can be fed into a Kalman 
filter in place measurement data from a gyroscope. The filter receives attitude update 
information from attitude sensors to produce smoothed attitude reference and angular rate 
error estimates to enhance satellite attitude determination accuracy. The plant error 
introduced into the filter is a combination of all the internal sensor and modeling errors 
that are used as inputs to the dynamic model. All known sensor biases must be 
incorporated into the plant model since the Kalman filter assumes errors to be zero mean 
Gaussian. Depending on different modes of operation, different plant noise models may 
be required for acceptable filter performance. Normally, filters that receive 
measurements from gyroscopes are designed to track gyro biases which remain relatively 
constant in the body frame. The rate error from the dynamic estimate, however, exhibits 
different characteristics. Since the dynamic model is external to the Kalman filter higher 
order effects can be more easily incorporated into the dynamics but onboard processing 


capability may still limit the complexity. 


I. PURPOSE 


The objective of this study is to develop and evaluate a practical attitude and rate 
estimation scheme for a multi-body spacecraft attitude control system that incorporates 
both real time angular rate calculation from the system dynamic model (dynamic gyro) 
and a Kalman filter estimator with attitude sensor updates provided by star trackers. It is 
hypothesized that acceptable attitude control performance can be realized by 
maneuvering multi-body spacecraft without hardware gyroscope data using this 


methodology. 


The evaluation of this approach is performed through simulation using a model 
developed in SIMULINK. A MATLAB script file is used to set up the necessary 
initialization parameters and system specifications. Simulation results are presented 


graphically in MATLAB plots. The performance of the developed gyroless attitude 
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determination system is compared to a conventional gyro-based system that uses the 
same Kalman filter estimation algorithm and attitude updates. Results are also presented 
that analyze the affects of several major error sources on the performance of this dynamic 


gyro based attitude determination system. 
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Hl. SIMULATION OVERVIEW 


The evaluation and analysis of attitude determination scheme proposed in this 
thesis is done through modeling and simulation. The context for the development of the 
attitude determination system is a multi-body attitude control system for a maneuvering 
spacecraft. An attitude simulation model for the spacecraft that includes vehicle 
dynamics, determination and control subsystems as well as modeled error sources is 
developed using MATLAB/SIMULINK. An overview of the simulation and top level 
subsystems is presented in this chapter. Subsequent chapters provide the descriptions of 
the subsystems and derivation of the equations on which the model is based. In the 
following chapters, actual SIMULINK diagrams are presented that show the 


implementation of associated subsystems. 


A. CONCEPTUAL ATTITUDE SIMULATION MODEL 


A conceptual representation of the simulation developed for this study is 


illustrated in Figure 2. 
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The model is divided into subsystems represented in color shaded blocks. This 
breakdown reduces the complexity of the overall model into manageable segments to aid 
in design and analysis. Arrows indicate the dynamic coupling and flow of data between 
subsystem blocks. The top level SIMULINK diagram that implements the concept is 


shown in Figure 3. 
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Figure 3. Top Level SIMULINK Attitude Simulation Model 


B. SPACECRAFT MOTION 


The actual spacecraft attitude motion is simulated in the rotational dynamics and 
kinematics subsystems with inputs and outputs represented by solid black lines. Relative 
motion of the secondary body or appendage is treated in a separate subsystem with 
dynamic effects coupled into the overall spacecraft motion through momentum exchange 
directly related to the drive motor rate. Dynamic effects of reaction wheel control are 


also realized through momentum exchange. Magnetic control effects are input to the 
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spacecraft dynamics as external torques along with modeled and approximated 
disturbance torques. The relative appendage motion also causes changes in the 
spacecraft’s moments of inertia used in the dynamic model. The kinematics subsystem 
propagates the actual spacecraft attitude quaternions referenced to an inertially fixed 


coordinate frame. 


C. CONTROL SYSTEM 


Primary spacecraft attitude control is conducted through momentum exchange 
with reaction wheels. The reaction wheel commands are based on control laws and up- 
linked feed forward torque command profile. The control laws use errors between the 
measured and commanded (desired) spacecraft angular rates and attitude quaternions. 
The relative momentum generated in the wheels is subtracted from the spacecraft 
momentum in the dynamics subsystem. Magnetic torques are generated based on 
reaction wheel momentum build up when the system is set for momentum dumping. 


These torques are directly applied in the spacecraft dynamics block. 


D. ATTITUDE DETERMINATION SYSTEM 


The attitude determination subsystem is of primary concern in this thesis. It 
incorporates the dynamic modeling concept and an error state Kalman filter in order to 
correct for attitude propagation errors. The option of using modeled mechanical gyros to 
determine spacecraft angular rate is also maintained in order to allow comparison 
simulations to be conducted. The concept of using a dynamic model for analytical rate 


determination will be referred to in the sequel as the dynamic gyro. 


The attitude determination computer simulation uses a trigger to control the 
bandwidth of discrete calculations. The basic data flow within the attitude determination 
computer is shown in the large rectangular subsystem block in Figure 2. The equivalent 


SIMULINK diagram that implements this subsystem is shown in Figure 4. 
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Figure 4. SIMULINK Attitude Determination Subsystem 


Random noise and drift rates are added to the actual angular velocity vector for 
the simulated rate input of the mechanical gyro option. Bias corrections from the Kalman 
estimator are applied to the angular rate before it is used in the attitude propagator. When 
the dynamic gyro option is simulated the angular rate is analytically determined from 
measured and known spacecraft parameters fed into a discrete dynamic model. Artificial 
errors and noise are applied to measured and derived parameters though internal sensor 
models. Momentum corrections to the dynamic gyro are derived from Kalman filter 


updates. The calculated rate is then supplied to the attitude propagator. 


The error state Kalman filter algorithm depends on updates generated by star 
tracker measurements. A star tracker model is used to produce artificial noise corrupted 
horizontal and vertical measurements related to the star tracker orientation. The Kalman 
filter tracks rate errors for bias correction of gyro measurements or momentum 
corrections for the dynamic gyro. It also produces attitude corrections that are applied to 


the attitude propagator output. 
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The attitude Propagator uses discrete kinematics to convert angular rate to 
quaternion attitude. These estimated parameter are then used by the attitude control 


algorithm to complete the feedback loop. 
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Hl. DYNAMICS AND CONTROL 


As three-axis stabilized spacecraft become more technologically advanced, their 
operations require more slewing maneuvers and their dynamic complexities increase. 
This often leads to multiple rigid or flexible components that have independent pointing 
and tracking requirements. Complex satellites often consist of a primary payload that 
demands strict pointing control while directional telemetry and command antennas or 
secondary payloads are controlled independently for other functions such as tracking a 


ground station throughout its maneuvers. 


Here we consider the dynamics of a spacecraft that consists of a primary body 
with momentum exchange control devices and a coupled rigid secondary body or 
appendage. The secondary body rotates with one degree of freedom relative to the 
spacecraft about an axis through the centers of mass of both bodies. Under these 
conditions, the position of the spacecraft center of mass remains stationary during 


appendage relative motion. 


A. BIFOCAL RELAY MIRROR SATELLITE DYNAMICS 


The example spacecraft for which the equations of motion are derived is the 
Bifocal Relay Mirror satellite shown in Figure 5. This spacecraft is designed to receive 
laser energy from a ground station through a receive telescope and to redirect it through 
an optically coupled transmit telescope to a different point on the ground. The primary 
body is the transmit telescope and the rigidly attached spacecraft bus while the smaller 
receive telescope is treated as the secondary body or appendage. Weight and balance 
design ensures that the centers of mass of the two bodies lie close to the coupled axis of 
rotation so that the system center of mass is nearly constant during relative motion. 
Pointing control of the spacecraft is accomplished with reaction wheels while a drive 


motor is used to control the relative angle between the transmit and receive telescopes. 
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Figure 5. Bifocal Relay Mirror Satellite 
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There are three coordinate systems defined for the development of the dynamic 
equations. These coordinate frames are depicted in the system diagram shown in figure 
x. The xyz coordinate system is fixed to the primary body with the origin at its center of 
mass. The x-axis is oriented parallel to the axis of rotation between the primary and 
secondary bodies and the z axis is parallel to the optical axis of the telescope. The y axis 
is defined such that the xyz coordinate system is a right-handed mutually orthogonal 
frame. The x’y’z’ coordinate system is similarly oriented to the secondary body with its 
origin at the center of mass. The x and x’ axes remain parallel during appendage motion 
while the angular displacements of the y’ and z’ axes from the y and z axes respectively 
are defined by the relative rotation angle 4. The equations of motion are derived for the 
spacecraft body coordinate system, XYZ, which is parallel to xyz frame. Its origin is at 


the total spacecraft center of mass. Unit vectors along the X, Y, and Z spacecraft body 


axes are given by 7, j and k respectively. 


L, Rigid Body Equations of Motion 


The angular equations of motion are derived from the application of Newton’s 
second law to rotational dynamics. In the general case the equation of motion defined in 


an inertial frame for a rigid body about an arbitrary point P is given by 
M, =H, + f(* «dn (3.1) 


where M, is the net external torque applied to the body about P, H, is the angular 
momentum of the body about P, fr, is the position of the center of mass relative to P, fi, is 


the vector from the center of mass to the position of dm in the body, and dm is an 


incremental unit of mass within the body 


If point P is coincident with the body center of mass, then fi, = 0 and the equation 


simplifies to 


M=H (3.2) 
Equation (3.2) applies in inertial reference frames. To extend it to body coordinates 


where it can be employed, we need to understand the concept of vector derivatives in a 
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rotating coordinate frame. If the body frame rotates relative to an inertial reference frame 


= 


with angular velocity U then the derivative of any vector A in inertial coordinates can 


be related to the derivative in the rotating body coordinates by 








+i xAl (3.3) 


Applying this relation to the angular momentum, the equation of motion for a rigid body 


in rotating body coordinates becomes 
M=H+ ix H G.4) 


De Multi-Body Equations of Motion 


For the two body Bifocal Relay Mirror system with reaction wheels shown in 


figure 1, the total system angular momentum in body XYZ coordinates can be written 


H=H, +H, +H, (3.5) 


where H_ is the total system angular momentum of the spacecraft, H, is the angular 
momentum of the system due to the rotation of the body coordinate frame, neglecting the 
contribution due to relative motion of the receive telescope and reaction wheels, H,,, is 
the angular momentum due to the relative motion of the secondary body, and H,, is the 
angular momentum due to the relative motion of the reaction wheels 

Since the Bifocal Relay Mirror spacecraft is approximated as a system of rigid 


bodies, we can substitute the total system angular momentum, H, into Equation (3.4) to 


get the multi-body spacecraft equation of motion 


M =H, +H,,+H, +0x(H, + H+ H,) (3.6) 


where M_ is the net external torque applied to the spacecraft about its center of mass 


including all control and disturbance torques. 
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3. Moments of Inertia 


To determine these angular momentums we need to consider the moments of 
inertia of the spacecraft and it components. We define I, to be the inertia matrix for the 


primary body (transmit telescope and bus) about its center of mass in the xyz coordinate 


frame. 
T T T 
Ix aly a 
L=|-l, ly Ty (3.7) 
T T T 
ag oli ie 


The vector from the center of mass of the primary body to the spacecraft center of mass is 


given by 
=X, +y,j+Z,k (3.8) 


The moment of inertia matrix for the receive telescope, k, about its center of mass given 


in the x’y’z’ coordinate frame is given by 


lbs Try TS 

_| qR R R 
Ip=|-FR, AR, (3.9) 

Ty, Ey leg 


The vector from the center of mass of the receive telescope to the spacecraft center of 


mass is given by 
Tp =X +¥ gi tZgk (3.10) 


To align the receive telescope inertia matrix with the spacecraft body coordinate 
frame a transformation matrix is applied. This time varying matrix is the direction cosine 
matrix that defines a single rotation. The axis of rotation is parallel to the body X axis 


and is of magnitude 4. The transformation matrix is given by 


1 0 0 
T, =|0 cos(a)  sin(a) (3.11) 
QO -sin(4) cos(a) 
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The SIMULINK subsystem used to generate this x-axis rotation matrix given an input 


angle is shown in Figure 6. 








Figure 6. SIMULINK Subsystem Diagram: X-Axis Rotation Matrix 


To obtain the total system inertia matrix about the spacecraft center of mass in the 
body coordinate frame XYZ the parallel axis theorem is applied to the inertia matricies of 
each body. Additionally the rotation matrix is applied to the receive telescope inertia 
matrix to align it with the body coordinates. If mp and mp are the masses of the transmit 


telescope and receive telescope respectively the total system inertia matrix is given by 


2,.2 
YrtZp -XqYp TX qZy 
= 2 2 T 
f= Let) Hxey- Zee. Hee |e Le 

2 2 

TXpZp 0 TY_qZp XP + Yt 
fw i (3.12) 

YrtZae —XrYr —XrZR 

2 2 

+Mpg|]—XpYp ZetXp —YrZpR 

2 2 

—XpZp —YrZp XptYr 


The calculation of the spacecraft system inertia matrix is implemented in the SIMULINK 


subsystem shown in Figure 7. 
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MOMENT OF INERTIA CALCULATION FROM APPENDAGE ANGLE 











Figure 7. SIMULINK Subsystem Diagram: Spacecraft Moment of Inertia Matrix 
a. Rate of Change of Spacecraft Inertia Matrix 


All components of the system moment of inertia matrix are constant except 
Ty1,T, which varies with the relative rotation angle 4. Therefore the rate of change of 


the system inertia matrix is given by 


I= Tl1,T, + TrlyTy (3.13) 
0 0 0 0 0 0 
where T, =|0 -sin(4)A cos(4)A | and Ty =]0 -sin(a)d_ -cos(4)4 |. It can be 
Q -cos(4)a  —sin(a)a QO cos(4)a  —sin(a)a 


shown that I as a function of 4 and 4 is given by 


0 
l@ 1] Sin(4) Bb €ps(a) 
Ix, sin(4) — Ky cos(4) 
Ixy sin(4)+ Ik, cos(a) 
21%, sin(4)cos(4) + 21%, cos*(4)— 21}, sin*(4)— 21}, sin(4)cos(4) 
Thy sin’( 4} Y, cos (4+ 4], sin(4)co@i)+1},cos’(4)— Tj, sin’(4) 
Is, sin(4) — Ih, cos(4) 
Thy sin’( 4} Y, cos (4+ 4], sin(4)cos(4} J, cos (4} J, sih 4) 
21%, sin’( 4+ 2, sin(4)cos(4} 2}, sin(a)cos(4} 2}, cs’( 4) 
(3.14) 


The SIMULINK subsystem diagram that calculates the rate of change of the spacecraft 
inertia matrix based on appendage relative angular orientation and rate is shown in Figure 
8. 
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COMPUTE DERIVATIVE OF SPACECRAFT INERTIA 
MATRIX FROM APPENDAGE ANGLE AND RATE 


column 
vector 








Figure 8. SIMULINK Subsystem Diagram: Rate of Change of Spacecraft Inertia 
Matrix 


4. Angular Velocity 


The spacecraft angular velocity is defined by the angular velocity of the primary 


body with respect to the inertial reference frame. This angular velocity vector is given by 
u,; =| (3.15) 


The relative angular velocity between the primary and secondary bodies due to the 


rotation about the body x axis is given by 


U rel 


a 
=| 0 (3.16) 
0 
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The angular velocity of the receive telescope is equal to the spacecraft angular velocity 


plus the relative velocity 


u,t+a 
Up =u; tt, =| uy (3.17) 
u, 
5; Angular Momentum 


The components of the total spacecraft angular momentums given in Equation 
(3.5) can now be defined. The system angular momentum neglecting the momentum due 


to the relative motion of the receive telescope and the reaction wheels is given by 
H, =I" (3.18) 
The angular momentum of the receive telescope relative to the spacecraft is given by 


A. =le wil (3.19) 


rel 


Substituting Equations (3.18) and (3.19) into Equation (3.5) we get the total system 


angular momentum 
H=I-4#,1j4 J (3.20) 


Substituting this relation into Equation (3.6), the spacecraft equation of motion can be 


rewritten as 
M=I--8 7 Wp Ig+) HO “Tg (3.21) 


6. Solving for Spacecraft Angular Rate 
To solve for the spacecraft angular rate, i, we can isolate i in Equation (3.21) 
and perform the integration using a computer solver. This however places an 


unnecessary burden on the processor to continuously calculate the time derivative of the 
spacecraft moment of inertia. A simpler method is accomplished by first solving for Hy . 


The spacecraft equation of motion is rewritten as 
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M=H, +1, Ha (3.22) 
which can be solved for H, ‘ 


H, =M-I, y- H+ Kw (3.23) 
After the integration, 
Ay = [(M-I, Q- H+ Fa) de f( Me Ha) de [W-H (3.24) 
t 


t 


the angular rate, U , is obtained from 
i =I'H, (3.25) 
Figure 9 shows the SIMULINK diagram for implements the spacecraft dynamics for the 


Bifocal Relay Mirror attitude simulation. The integration is performed using the 


Dormand-Prince ode5 solver. 


CONTINUOUS 
SPACECRAFT 
DYNAMICS 





Figure 9. SIMULINK Subsystem Diagram: Spacecraft Dynamics 
B. COMMAND 


Maneuvering spacecraft often require feed forward command to maintain precise 
tracking requirements throughout their maneuvers. The Bifocal Relay Mirror satellite 


must maintain tightly controlled ground tracking by both the transmit and receive 
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telescopes in order to perform its mission. For this system the envisioned feed forward 
command will include control torques necessary to maintain a dynamic attitude profile, 
which includes the relative motion between the primary and secondary bodies, calculated 
in the absence of disturbances. The spacecraft control laws will be used to correct for 
errors in the calculated command profile and external disturbances based on errors in 


spacecraft attitude and angular rate from the attitude determination system. 


The estimated spacecraft rotational profile is predetermined prior to the execution 
of a maneuver and tracking operation. The profile includes the spacecraft body attitude, 
angular rate and angular acceleration as well as the relative angle, rate and acceleration 
between the transmit and receive telescopes. The net external torque required to maintain 
this profile in the absence of disturbances is fed forward to the control devices. This 
torque can be calculated from the spacecraft equation of motion, Equation(3.21), 


neglecting the reaction wheel control devices 
M.=1°, 8 Bg hdt ot 7B Ted) (3.26) 
In this equation the subscript c is used to represent feed forward command. The 


equivalent SIMULINK subsystem used in the simulation to generate the feed forward 


command torque is shown in Figure 10. 
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Figure 10. SIMULINK Subsystem Diagram: Feed Forward Torque Command 
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C. CONTROL 


Primary attitude control of the Bifocal Relay Mirror satellite is accomplished with 
reaction wheels. Four reaction wheels are arranged in a pyramid constellation to achieve 
redundancy in the event of a single wheel failure. Under normal operations three wheels 
are operating while the forth is off line. Control torques commands are calculated in 


body coordinates and then distributed to the three operating wheels. 


Torque commands are generated from the feed forward command profile plus the 
attitude control laws with compensation for the gyroscopic torques generated from the 
spinning reaction wheels. The simple partial plus derivative type controller is chosen for 
this satellite. The control laws are based on attitude quaternion and angular rate error as 
calculated from the outputs of the attitude determination system where the quaternion 


error is calculated by Equation (4.16) [Ref. 9]. The three body axis wheel control laws 


are given by 
Hes = K det Ki u 
Hy = K dent K, uU (3.27) 
Hy3= K dest K, u 


More exotic control schemes exist but the optimal state estimates provided by the 


attitude determination system are simply employed by these control laws. 


The control law gains are chosen to minimize steady state errors attitude errors 
while ensuring that oscillations induced by the attitude control system do not interfere 
with on-orbit structural modes and payload components. To determine the optimal 
control gains, attitude control simulations are conducted with a_ representative 
maneuvering profiles and external disturbances. The state inputs used by the controller 
during gain adjustment simulations are perfect attitude and rate knowledge. As explained 
in Chapter I, the optimal controller can be determined independently from the attitude 


determination system since its outputs are based on optimal state estimation. 


The reaction wheel control subsystem used in the SIMULINK attitude simulation 
model is shown in Figure 11. Saturation and rate limiting is applied to simulate nonlinear 


effects in real reaction wheel control systems. The net control torque is determined and 
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then distributed to the individual reaction wheels based on their orientation within the 


constellation. 


REACTION WHEEL CONTROL 


-Kq1*(u[14)}* u[4}+u[13) ufS)-u[12}* u [6)-u [1 1)" u [7] Ke *Cu[1]-ufep 
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Figure 11. SIMULINK Subsystem Diagram: Reaction Wheel Control 
D. DISTURBANCE TORQUES 


Representative disturbance torques are simulated in order to observe spacecraft 


attitude performance in a realistic space environment. 


1. Gravity Gradient Torque 


At low altitudes the torque induced by gravity gradient on spacecraft without 
matched body moments of inertia can be significant. A model for the gravity gradient 


torque on a spacecraft in orbit is given by 


Be 3i - = 

M, Spe dere (3.28) 
0 

where Ry the distance from the spacecraft center of mass to the Earth’s center, R, is the 


unit vector in that direction given in body coordinates, and I is the total spacecraft 
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moment of inertia matrix. If the direction cosine matrix from the orbit to body frame, 


®C° , is known, R, is given by 


0 C3 
Ry = ?C°| 0} =| 43 (3.29) 
1 C45 


In the attitude simulation model, the orbital reference frame direction cosine matrix 
components are propagated using the spacecraft angular rate relative to the rotating orbit 
frame from an initial orientation. The relative orbital rate is determined from the inertial 
rate minus the rate of rotation of the orbit frame. The SIMULINK subsystem for 


propagating the orbital reference is shown in Figure 12. 





qo ORBITAL REFERENCE CALCULATIONS 








Figure 12. SIMULINK Subsystem Diagram: Orbital Reference Propagation 


Using the orbital reference coordinates the gravity gradient torque can then be 


written 
31 C13 Ci3 
0 
C33 C33 


Figure 13 shows the SIMULINK subsystem that models the gravity gradient torque based 


on the orbital reference frame direction cosine matrix. 





GRAVITY GRADIENT CALCULATION FROM ORBITAL REFERENCE 











Figure 13. SIMULINK Subsystem Diagram: Gravity Gradient Torque Model 
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Ze Other External Disturbance Torques 


Other disturbance torques include those due to unbalanced solar pressure, 
magnetic interactions, and aerodynamic drag effects. These disturbances are not as easily 
modeled as gravity gradient but their characteristics and magnitudes are important 
considerations in the design of spacecraft attitude control systems. In the simulation 
model, secular and slowly varying periodic moments are introduced in each body axis to 
account for these unknown disturbances. Worst case magnitudes are chosen based on 


orbit profile and spacecraft characteristics to ensure robust control design. 


E. MAGNETIC MOMENTS 


The magnetic moment imparted on the spacecraft by the earth’s magnetic field is 


dependent upon the magnetic field strength, B, and the spacecraft’s magnetic dipole 


vector, m. The magnetic moment is given by 


M,, =mxB (3.31) 


m 


1. Magnetic Field Model 


Highly accurate models of the earth’s magnetic field have been developed but a 
simple dipole model is adequate for the purposes of this simulation. This approximation 
assumes a simple dipole magnetic field tilted 11 degrees from the earth’s spin axis. The 
earth’s magnetic field is a function of the earth’s unit dipole vector, M , the distance from 
the center of the earth to the center of mass of the spacecraft, R, and the unit vector in 


that direction, R. The magnetic field is given by 














B=—.[3(W R)R-M | (3.32) 


where K is the earth’s magnetic field constant equal to 7.943 x 10" Nm’/a’. Using 
classical orbital elements it can be shown that the components of the magnetic field 


vector in orbital coordinates are given by 
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cos )(cos(asin(i) —sin(a)cos( } cos( H )-sin( sin (a)sin( u) 
BB == —cos(a)cos(i) —sin(asin( } cos( 2 
2sin(v) (cos(sin() -sin(&cos( ¥cos( 9 )-2cos( ys in(dsin(w) 

(3.33) 
where v is the true anomaly of the spacecraft, € is the magnetic dipole tilt angle, 7 is the 
orbit inclination, and uw is the angle of the magnetic dipole normal to the line of ascending 
nodes. The SIMULINK subsystem that propagates the earth magnetic field vector with 


the orbital reference is shown in Figure 14. 


EARTH MAGNETIC 
FIELD CALCULATION 


parmtrs 


fe Bm*(coxtul[1]F(costul2]Fsin(ulS]}sin(ul2Z]F costu[s]F coxu[4}-sin(ult]Fsintul2}Fsin(u[4]}) i 
in Bm*C-costul2]F cox(u[3]}sin(u[2]Fsin(u[s]F coxtu[4]) bl 
Bm*(2*sin(u[t]P(comul2]Fsin(u[s]}}sin(u[2]F coxuls]F costul[4] +2" cos(u[1] Fsin(ul2]Fsin(u[4]) 








Figure 14. SIMULINK Subsystem Diagram: Earth Magnetic Field Vector in Orbital 
Coordinates 


The magnetic field vector can be transformed to the spacecraft body coordinates 


with the direction cosine matrix. 
B, = °C°B, (3.34) 


where °C° is the transformation matrix from the orbit coordinate frame to the body 
frame. Figure 15 shows the subsystem that realizes this transformation using direction 


cosine matrix components. 
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Figure 15. SIMULINK Subsystem Diagram: Magnetic Field Conversion to Body 
Coordinates 


2. Magnetic Control Torque 


Using magnetic torquers a spacecraft magnetic dipole can be generated to react 
with the earth’s magnetic field to produce a control torque. In this simulation magnetic 
control can be used to help desaturate the reaction wheels. The control law for magnetic 


dumping of reaction wheel momentum is given by 


Mona =—k(B, XH, ) (3.35) 


where k is the magnetic torquer gain. Figure 16 shows the SIMULINK subsystem that 
simulates the magnetic dumping control torque. Saturation is added to simulate nonlinear 


effects in the torque rods. 
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Figure 16. SIMULINK Subsystem Diagram: Magnetic Torquers for Momentum 
Dumping 


Magnetic torquers can also be used for attitude control. The control laws for the 


torquers to produce a desired control moment, M,,,, on the spacecraft is given by 


MoS er (3.36) 
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IV. ATTITUDE REPRESENTATION AND KINEMATICS 


There are many ways of describing the orientation of one coordinate system 
relative to another. The most common descriptors used in spacecraft attitude 
determination include Euler angles, direction cosine matricies and quaternions, also 
known as Euler parameters [Ref. 10]. Euler angles provide a convenient way to represent 
attitude and are usually he easiest to visualize. However, singularities arise when the 
relative orientation from the reference coordinate system becomes large. Therefore, 
highly maneuverable spacecraft require other means of attitude representation. Direction 
cosine matricies and quaternions overcome this problem. Direction cosine matricies 
provide the most convenient way of transforming vectors between coordinate systems but 
require significantly higher attitude processor throughputs than quaternions. As an 
example, propagating an attitude matrix with angular rate data requires the integration of 
nine elements while the quaternion has only four. For these reasons, the quaternion 


representation is chosen for this attitude simulation model. 


A. QUATERNION DEFINITION AND CHARACTERISTICS 


The four vector quaternion based representation G¢ R* can be written as 


q= (4.1) 


or equivalently, 
q=qitqy tak +q, (4.2) 
where the quaternion has a three vector imaginary part 


qi 
q; =| q, |and a scalar real part qy =q,. 


q3 
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1, Meaning of Quaterions 
Quaternions represent the angular orientation of a body relative to a reference 
coordinate frame by a single axis rotation of magnitude é about the eigenvector axis 


given by é corresponding to the +1 eigenvalue of the direction cosine or attitude matrix. 


The direction cosine matrix can be written in terms of these four parameters 


10 0 0. = ge 
A=|0 1 O|cos(®)+| é; 0  -€, |sin(@)+&™(1 —cos(@)) (4.3) (4.4) 
001 ce. § “e| er 


Similarly these parameters can be represented by the four quaternions 


és sles ee 
Gi =| 942 |=P 291N (=| 
gg 3 (4.5) 


= =COS Z 
Gr =4 5 


which have the added property that 


ja] =a5 +43 +95 +4; =1 (4.6) 


2. Attitude Matrix 


The equivalent attitude or direction cosine matrix can be generated from 


quaternions using 


G=q—=a+q, 2Gq+aq) -2(a4,—0,4;) 
A=| 2(q,q.-4394) a; +q3-q3t+q;  2(q.q3+49)94) (4.7) 
2(qiq3+4.04)  2(q.q;-4)q,) Gq; -G3 +aq3 +4; 


The SIMULINK subsystem used in this simulation model to convert a quaternion to the 


equivalent attitude matrix is shown in Figure 17. 
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Figure 17. SIMULINK Subsystem Diagram: Attitude Matrix from Quaternions 


3: Quaternion Multiplication 

The multiplication of two quaternions to form a third defines an angular 
orientation resulting from two eigen-axis rotations. If q represents the transformation 
from coordinate frame A to B and q’ represents the transformation from frame B to C 
then the transformation from frame A to C is given by 

q=qq° (4.8) 

This multiplication can be implemented several ways. The two quaternions can 


be multiplied directly using the quaternion format given in Equation (4.2) 
(qita,j +q3k+q,)=(qita, i+ ak+ qi qt q3j+ ak +4) (4.9) 


using the equalities 


P=j*=k?=-1 

ie ae (4.10) 
jk=—-K =i 

ki=—ik = j 


Quaternion multiplication can also be performed by treating the imaginary and real parts 


separately 
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qi= qeai+ ard yl xqr (4.11) 
Gr = rr ~ 441 





In matrix form quaternion multiplication is given by 


aq} }ds -3 4 {Iq 

42/_]% d+ ha] a (4.12) 
q3 “G2 di; Ga 3 |} G3 
dal [di 2 43 a]lqa 


This multiplication method is used in the attitude simulation. The SIMULINK subsystem 


is shown in Figure 18. 





QUATERNION MULTIPLICATION g=q'q" 


fe u[4fufS)-usfuf6}+ul2fufF)}+u[1} us] ei 
uf3fu5)+4 [47 u[6)-ult ul? )}+u[2 us] 
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Figure 18. SIMULINK Subsystem Diagram: Quaternion Multiplication 


4. Coordinate Rotations 


Sequential axis rotations such as those that define the three Euler angle 
representation can be realized by the successive quaternion products. Additionally, if 


there are three small simultaneous rotations é,, €,, and €, about the coordinate axes x, y, 


and z respectively the resulting quaternion is determined as follows: 


eC, 
Let @= (8; +8; +83)'” and 6-8 , |, then 
C3 
qi 
gq, =| 4 =6.in( 2) and qz =a,=00($ (4.13) 


a. Quaternion Inverse and Identity 


The inverse of a quaternion is its complex conjugate and is analogous to the 
transpose of a direction cosine matrix. A quaternion is conjugated by reversing the sign 


on the vector part. 


q, -q\ 
li-q= do then q = “da 
q3 “3 
q4 q4 


The identity quaternion is found by multiplying any quaternion by its conjugate as 


shown below. 


(4.14) 


6. Quaternion Error 


A difference quaternion can be defined between to orientations that are referenced 
to the same coordinate frame. If q’ represents the transformation from coordinate frame 
A to B and q’ represents the transformation from frame A to C then the transformation 


from frame B to C is given by 


q=qqg =| , i < a eee (4.15) 


This can be used to calculate the error quaternion where the target orientation is 


given by q and the actual spacecraft orientation is q’ . 


=8 pki 


q-=94 4 (4.16) 
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The SIMULINK subsystem diagram for determining the error quaternion is given in 


Figure 19. 


QUATERNION ERROR ge=q'*q" 


fe ufafufS}+ufsPuf6)}-uf2Pul?)}-u[t Pus) al 
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Figure 19. SIMULINK Subsystem Diagram: Error Quaternion 


The angular difference in radians between the original quaternions is contained in 


the real part of the error quaternion. 


&, =2cos |(q,.) (4.17) 


7. Vector Transformations 


Transforming vectors between coordinate systems requires two quaternion 
multiplications where the vector is treated as quaternion with a real part of zero. If q 
represents the orientation of coordinate frame B with respect to reference frame A then a 
vector given in coordinates of the reference frame, V,, can be transformed to the 


coordinates of frame B by 
Vp =q VAG (4.18) 


B. QUATERNION KINEMATICS 


Actual spacecraft motion is simulated using the continuous quaternion differential 
equations. In the attitude determination computer simulation the attitude propagator uses 


the discrete kinematic equations. 
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1. Continuous Kinematics 


The differential equation for the quaternions of a rotating coordinate system can 
be found by differentiating this equation with respect to a fixed reference frame A. If the 


rotational rates for coordinate frame B are given by 














uy 
dv dv dv ho 
ii =| a, |,then —4) =—#4) =0 and —#| =uxV, 
2 dt |, dt} I 
u 








The inertial derivative of equation (4.18) becomes 


(4.19) 


Using the properties of quaternions, substituting the original equation for Vv, and 


noting that the vector V, is arbitrary it can be shown that 


qii (4.20) 


qi qi Ou a) 11} 
jae = = =) O ss. NX 
dq _ 4a wl qo nas ald 2 u; = qo (4.21) 
dt q3 2 qs 2}u, a, O U,]I/q, 
a qa -u, -U, -u;, 0 |i q, 
= S(t) a eae ; : ; : 
where S(U)=| _- 0 and S(U)is the skew symmetric matrix associated with the 
-u 


vector U. The kinematics subsystem diagram used in the SIMULINK model is shown in 


Figure 20. The integration is performed using the Dormand-Prince ode5 solver. 
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Figure 20. SIMULINK Subsystem Diagram: Continuous Kinematics 
Zs Discrete Kinematics 


The solution to the continuous differential equation when U is held constant is 





1st _ ine 
q= AC hi = cosé \Tiga +S Es ne ) “0 (4.22) 


26 0) 


N 


where €= 











This solution leads to the discrete implementation of the quaternion 


kinematic equations. 

The attitude determination system uses angular rate information to kinematically 
propagate the spacecraft attitude quaternion in discrete time steps. The angular rate 
vector ie R° is integrated with time step At to produce the incremental angle vector 


Aéeé R®. For small time steps Aé approximates an eigen-axis rotation in the current 


body coordinate frame so it can be related to the change in the attitude quaternion by 








“ ei ies seo EN (4.23) 


The SIMULINK subsystem used in the discrete attitude propagator to obtain the 


incremental quaternion step is shown in Figure 21. 
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Figure 21. SIMULINK Subsystem Diagram: Discrete Attitude Propagator 


The new attitude quaternion is simply determined by 


Ges =o ( A (4.24) 


where the equation implies quaternion multiplication. 
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Vv. KALMAN FILTER 


The Kalman filter provides a non-deterministic means of estimating a system state 
vector using a state space mathematical plant model and sensor measurement data related 
to some subset of the state variables. It is a stochastic optimal estimator designed to 
minimize the weighted mean square error in the state estimate. The Kalman gain matrix 
determines the weighting based on the relative confidence between the past state estimate 
propagated to the current time and the current partial measurement of state variables. 
The error covariance matrix, the second statistical moment of the state vector, tracks the 
confidence in the state estimate while a measurement error covariance matrix relates the 


confidence in the measurement. 


A. RECURSIVE DISCRETE KALMAN FILTER 


The Kalman filter can be implemented recursively since all of the information 
from past state measurements is encapsulated in the previous state estimate and error 
covariance matrix with are both tracked. During state propagation the error covariance is 
updated to reflect additional error added by imperfections in the plant model. A recursive 
discrete Kalman filter is used in the proposed attitude determination scheme. Figue 22 


illustrates the recursive nature of the Kalman filter. 
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Figure 22. Discrete Kalman Filter Loop 


The standard equations for a recursive discrete Kalman filter are summarized in 
Table 2 below. For full background development and equation derivations see Appendix 


A [Ref. 1]. 





System/plant Model XO=X Wt 


W, ~ NO,Q,) 
Measurement Model Z, =H,X, +V, 
v, ~NOR,) 
Initial Conditions x HEE 





P, =E[(¥p —Xy )(Kp—Xp) 1 
Assumptions E[w,v,]=0 
(uncorrelated errors) forall. & 


Prediction State Estimate Extrapolation XO=x,1 
kl kl 
Error Covariance Extrapolation PO=P,66, (Qt a 
Correction Kalman Gain Matrix K, = P-Hi (H,P,;H;, +R ay 
State Estimate Update xi =x_+K,(Z,-H,x;) 
Error Covariance Update Py = (I- K,H,)P, 
Table 2. Discrete Kalman Filter Equations 
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B. ERROR STATE EXTENDED KALMAN FILTER IMPLEMENTATION 


A linear error state extended discrete Kalman filter is implemented in the 
simulated attitude determination model to estimate spacecraft attitude and angular rate. 
The nonlinear attitude propagation is performed discretely outside of the Kalman filter 
according to Equations (4.23) and (4.24) at the high bandwidth frequency of the attitude 
processor. The angular rate estimate used in the kinematic model is provided by the 
spacecraft gyros or pseudo gyro rate calculations. The Kalman filter is designed to 
provide bias corrections to the gyro outputs for attitude propagation [Ref.11]. If the 
pseudo gyro rate is used, the bias error is treated as a rate error from which a correction to 
the system angular momentum is determined. Measurement updates are provided by star 
trackers. The difference in the measured and predicted star vector is related to the 


attitude error and used to provide corrections to the state estimate. 


1. State Variables 


The Kalman filter used in this model estimates six state variables: a three vector 


of attitude errors, &€ R°, and a three vector of gyro bias errors, be R*. The total state 


vector is given by 
a Feo ae 
x=|_|ER (5.1) 


The attitude error, é , represents the deviation in the spacecraft attitude relative to 


the inertial reference frame given by a vector of three simultaneous rotations. The bias 


error, b, represents the change in bias of the angular rate data in the spacecraft body 


coordinate frame. 


2s Attitude Propagation Error Correction Methods 


There are two ways to implement the Kalman filter corrections to the attitude 
propagator. In the first method, the attitude propagator is fed with raw gyro rate 


information and the Kalman filter maintains the total gyro bias tracking with the state 


bias error, b. In this method, the attitude quaternion must be corrected at each time step 
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by the filter. The attitude error vector, é, is converted to an incremental quaternion 
rotation, AG and applied to the propagated attitude quaternion as shown in Chapter IV 
with Equations (4.23) and (4.24) or, using the small angle approximation Aq is given by 


Agq=—i+—j+—k+4+1 5.2 
q 5 re 5 (5.2) 


Then by quaternion multiplication the propagated attitude quaternion is updated 
Cis =(AQ)Guia (5:3) 
The alternative method of applying the filter correction to the propagator is to use 


the bias error state, b, to correct the unknown gyro bias, U, each time a measurement 


update is obtained 


Upnew = Uo +0 (5.4) 
For this method, the total unknown gyro bias is tracked separately from the filter and 
added to the gyro angular rate before it is fed to the attitude propagator. 

A =(ii, +u,) At (5.5) 
The attitude error state correction is then only necessary at measurement updates. Using 
this method allows the Kalman filter error vector x to be reset after each bias 
measurement update. Since the filter approximates nonlinear errors with a linear model, 


keeping the errors as close to zero as possible improves the estimate. This Kalman filter 


correction method is used in the simulated attitude determination model. 


3. Plant Model 


A nominal plant model is chosen for this Kalman filter which assumes a constant 
rate bias between filter updates and uncoupled (linear) attitude errors between 


propagation steps. It is derived from the continuous equation of state 


dx_|dé/dt} [o,, A” 
dt db/ dt 03x3  O3x3 |b 


where A is the estimate of the spacecraft attitude or direction cosine matrix. This linear 


| =f (t)x (5.6) 
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state equation includes coupled rotational effects through the time variant estimated 


attitude matrix which changes dynamically. 


4. State Transition Matrix 
For the discrete Kalman filter implementation of this model the state transition 
matrix, 6,, 1s calculated for each time step using the matrix exponential of f(t) from the 
continuous state equation 
Bost I A t 
6, = ef WAt -| 3X3 AA | (5.7) 
035 3 I 3X3 


where At is the transpose of the estimated attitude matrix at the current discrete time, tk 


and At is the interval between steps. The SIMULINK subsystem used to calculate the 
state transition matrix is shown in Figure 23. The estimated attitude matrix is calculated 


from the propagated attitude quaternion using Equation (4.7). 


STATE TRANSITION MATRIX 


Adt Trans 








Figure 23. SIMULINK Subsystem Diagram: State Transition Matrix 
5. Kalman Filter Prediction Equations 


The state transition matrix allows the state estimate, x, and associated error 
covariance matrix, P, to be propagated forward in the discrete Kalman filter prediction 


step 
Xin= Op (5.8) 


Pi= O,PO + (5.9) 
The subscript k+1 indicates next discrete time step, {,;, and the superscript — indicates 
predicted future value based only on information up to the current time step t. The 


covariance matrix relates the confidence in the associated state estimate. A larger P 
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indicates less confidence which means that the attitude determination system is assumed 
to have larger errors. Measurement updates taken when P is large will greater impact on 
the estimate than those taken when the confidence is high. The Kalman filter prediction 


step is implemented in the SIMULINK model as shown in Figure 24. 
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Figure 24. SIMULINK Subsystem Diagram: Kalman Filter Prediction Step 
6. Plant Noise Covariance 


The plant noise covariance, Q, is a positive definite matrix that characterizes the 
plant error accumulated through the time step, At, assuming that it can be modeled as 
Gaussian white noise meaning normally distributed with zero mean. Q is normally taken 
to be a diagonal matrix meaning that there is no known correlation between the errors of 
the six state variables. In terms of the variances of the state variables the plant noise 


covariance matrix is given by 


6a; 90, oe Oi. OY. 20 
06 30 0 0 O 
o- 0 0 4 O O O (5.10) 
0 0 0 ¢ O O 
0 0 0 0 4 O 
0 0 0 0 0 4 


If the time step is varied or the plant error is known to be significantly changed 
due to operating mode, then Q should be modified to reflect that change. Decreasing Q 
effectively asserts that the plant generates less error and the predicted state estimate will 
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be weighted higher relative to the measurement updates. On the other hand, if it is 
increased less confidence is placed on the estimate predicted by the model and more 
significant corrections will be applied during measurement updates. In this simulation, Q 


is taken to be constant. 


7. Kalman Filter Initialization 


In order to initiate the Kalman filter algorithm, an initial state estimate, x, and its 
associated covariance, P;, must be chosen. Least squares batch processing can be 
performed on two or more star tracker measurements prior to initiating the filter or an 
educated guess can be used. Since the state represents errors from truth, the selected 


initial state estimate is 


0 100 O O 0 0 0 
0 0 100 0 0 0 0 
-_|O) 0 0 100 0 0 0 
X,= with associated covariance P.=n . The 
0 0 0 0 1 0.0 
0 OO. Os Or OD 
0 0 0 0 00 1 


constant, n, is chosen to reflect confidence in the estimate at initialization. The attitude 
errors generated are normally about an order of magnitude up from the rate bias errors so 


the variances are weighted higher. 


8. Sensor Measurement Update 


The state estimate and the covariance matrix propagate forward without 
correction through each time step until a measurement update is produced. The state 
estimate is used to kinematically propagate the attitude quaternion while the covariance 


matrix builds up due to the added error through each step. 


a. Measurement Vector 


If an attitude measurement update is produced at the next time step, k+1, 
by one or more star trackers the Kalman filter correction step is applied. The star trackers 


produce horizontal and vertical outputs (H,V) corresponding to the position of the star on 
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the detector array. These outputs are internally compensated for with software and 
calibrated to form a measurement vector in tracker coordinates 
H —s 
a : S seietieak ae, 2S 
Sm=| V | with the normalization s,,=— (5.11) 


a 


which is the star tracker measurement vector. For simulation modeling purposes, noise is 





added to the H and V components to create an artificial star tracker measurement vector. 


The measurement vector is generated in the SIMULINK model as shown in Figure 25. 









ARTIFICIAUNOISY STAR 
MEASUREMENT WITH 
NORMALIZATION 





Figure 25. SIMULINK Subsystem Diagram: Star Tracker Measurement 


At every integration cycle, At, the attitude estimate is available from the 
Kalman filter prediction step. If a star tracker observation has occurred during that time 
step, the measurement vector must be propagated ahead to correspond to the current 
discrete time of the attitude computer. Otherwise, the state vector and the state transition 
and covariance matricies must be interpolated to accomplish the filter update. Tracker 
processing latencies and transport delays must also be compensated for in the time 
difference. For simplicity in this model, it is assumed that the all observations occur at an 
exact discrete time step of the integration cycle. To produce star measurements for the 
model, an artificial star tracker reference is chosen with representative noise applied to 


the H and V outputs. 


b. Predicted Measurement 


The predicted vector in star tracker coordinates is needed to determine the 
measurement residual for the update. This is generated using the known position of the 


star in inertial space. The detected star goes through the identification process and gets 
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compensated for aberration to yield a unit vector, g, in inertial coordinates. Using the 


estimate for the inertial to body attitude matrix generated from the filter state prediction, 


A 


A,,,» and the calibrated body to star tracker transformation matrix, T, a predicted vector 
is generated in tracker coordinates. 


A 


3=T A,,,3, (5.12) 


P 


The inertial star vectors for the model are generated by applying the 
transpose of the body to star tracker and true attitude matricies to the same star tracker 


reference used to generate the measurement. 


H 
S$=Aial | V (5.13) 
1 
Figure 26 illustrates the SIMULINK subsystem used to produce a simulated inertial star 


reference vector. 
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Figure 26. SIMULINK Subsystem Diagram: Inertial Star Vector 
c. Measurement Residual 


The measurement residual, z,;, 1s formed by subtracting the predicted 


from measured star vector and considering only the first two components. 
Zi41=E(S,, —S,) (5.14) 


0 


1 
where E= 
i 1 O 


| The measurement residual is determined in the SIMULINK model 


as shown in Figure 27. 
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Figure 27. SIMULINK Subsystem Diagram: Measurement Residual 
9. Observation Matrix 


Since the measurement residual is simply the difference between unit vectors, it is 
easily related to the state attitude errors. The observation or feedback sensitivity matrix 


defined by the relationship to the state variables 


Hyi=[ETA,.:$ (8) Oe (5.15) 


where S(S,)indicates the skew symmetric matrix associated with the inertial star 


reference vector. Figure 28 shows the subsystem diagram used in the SIMULINK model 


to produce the observation matrix. 
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Figure 28. SIMULINK Subsystem Diagram: Observation Matrix 
10. Kalman Gain 


The Kalman gain can then be calculated for the correction step at k+1 from the 


standard discrete filter equation 


= T - T -l 
Kyi =P A (Ryu +H Pe Hy. ) (5.16) 
where Re R*? is the measurement noise covariance matrix associated with the assumed 


Gaussian white noise in the H and V outputs of the star trackers. Since these errors are 


54 


usually similar and uncorrelated, the R matrix is normally taken to be a multiple of the 


identity. 


ae 6; 0 
06% 


Noisier or less accurate star trackers will generally have a larger R matrix which 
decreases the Kalman gain and thus provides less of a correction to the state and 
covariance matricies. In reality, the measurement error may not be constant for each 
correction step even when the same star tracker is used. Higher intensity stars, although 
easier to detect and identify, produce slightly larger noise. Also stars detected toward the 
field of view limits of the tracker usually have larger errors due to distortions than those 


detected near the star tracker optical axis. 


The Kalman gain SIMULINK subsystem is shown in Figure 29. 
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Figure 29. SIMULINK Subsystem Diagram: Kalman Gain 
11. Kalman Filter Correction Equations 


With the Kalman gain, the state and covariance matricies can be corrected with 


the measurement update 


A 


AXx=— Ky:Z,41 (5.17) 

Xj =XAKt (5.18) 
and 

Pri =C x6 mr Ky Ay) Prat (5.19) 


The covariance matrix update is accomplished in the simulation model as shown in 
Figure 30. 
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Figure 30. SIMULINK Subsystem Diagram: Covariance Updeate 


In general, the estimated state drifts away from the true state as it is propagated 
forward with the non-ideal plant model and is corrected back toward truth as 
measurement updates occur. Corrections to the attitude components associated with the 
body axes closely aligned with the star tracker optical axis will be corrected much less 
than those that are aligned perpendicular. The spacecraft attitude error continues to grow 
during periods where there are no cataloged stars in the sensor field of view or when the 
tracker information is unavailable. Additionally, when only one tracker provides updates 
for extended periods the angular orientation about that star tracker’s axis in body 


coordinates goes unchecked. 


Cc. CONTROLLER DESIGN IMPLICATIONS 


Since the Kalman filter is an optimal least squares estimator, the development of 
an optimal controller can be accomplished independently. Therefore, the attitude 
determination output should not affect the controller design. In this model, quaternion 
and rate error control laws are used as well as feed forward torque to generate reaction 
wheel commands. The optimal gains for the controller as determined by simulation with 
ideal deterministic attitude knowledge remain optimal with the attitude determination 


based on Kalman filter state estimates. 
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VI. DYNAMIC RATE CALCULATION 


The continuous dynamic equations of motion for the Bifocal Relay Mirror 
spacecraft are derived in Chapter III]. These equations produce the spacecraft angular rate 
from external control and disturbance moments applied to the body. A similar discrete 
model can be applied in the spacecraft attitude processor software to produce a real time 
calculated estimate of the angular rate, referred to as the dynamic gyro. This process is 
borrowed from The Aerospace Corporation’s “Pseudo Gyro” concept [Ref. 8]. At high 
bandwidth processor execution the discretization of the dynamics introduces little error. 
The angular rate generated by this method can be used as a substitution for conventional 
gyroscope outputs. Attitude determination based on the dynamic gyro can be 
implemented as a back up failure mode or a primary operating mode to increase the 


expected lifetime of the satellite gyroscopes. 
A. DISCRETE EQUATIONS OF MOTION 
The discretized equations of motion are derived from 
AH =)" MAt (6.1) 


where yM is the sum of external moments applied to the spacecraft including controls, 


modeled disturbances and gyroscopic stiffness. This allows the total system angular 


momentum to be tracked with 


H,,,, =H, + AH (6.2) 
Then subtracting the relative momentum of the reaction wheels and secondary body 
produces 

H, Hn (6.3) 


The calculated spacecraft angular rate is then given by 


i =r'A, (6.4) 
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B. MOMENTUM CORRECTION FROM KALMAN FILTER UPDATES 


The accuracy of the dynamic angular rate calculation ultimately depends on the 
tracking of the system angular momentum. If uncorrected, the numerous error sources in 
the model will cause the angular rate error to grow over time. The error state Kalman 


filter designed for gyro-based attitude determination systems can be used to provide the 


necessary model corrections. The gyro bias error states, b, are interpreted as spacecraft 
body rate errors. Using the calculated spacecraft inertia matrix, I, a correction to the 


system moment of inertia can be generated by 


AH... =Ib (6.5) 
The Kalman filter momentum correction is applied as if the error in the dynamic gyro is 
attributable to the total spacecraft body. The relative momentum terms from the 


secondary body and the reaction wheels are treated as if they are without error. 


The SIMULINK dynamic gyro subsystem used in the attitude simulation model is 
illustrated in Figure 31. 
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Figure 31. SIMULINK Subsystem Diagram: Dynamic Gyro 
C. INPUTS AND ERROR SOURCES 


After initial calibration, kinematic plant error in gyro-based attitude determination 
systems is almost entirely attributable to a single set of imperfect gyroscope rate sensors. 


As long as gyro data does not become erratic a Kalman estimator based on a slowly 
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changing rate bias plant model produces an effective attitude determination system even 


with relatively noisy rate inputs. 


The error in rate calculations from dynamic modeling, on the other hand, is due to 
numerous factors and is much harder to characterize. There are multiple internal sensors 
involved as well as dynamic modeling simplifications. Since the dynamic calculation is 
produced from total system momentum tracking any error in knowledge of external 
torques directly correlates to rate error. Errors in system or component moments of 
inertia have the same effect. Like gyro outputs, internal position and rate sensor data are 
corrupted by measurement and alignment errors. These data from all moving appendages 
and momentum exchange devices are critical to the accuracy of the rate calculation. 
Satellites not designed to use dynamic rate calculations for attitude determination are 
usually not equipped with appendage relative rate measurement sensors. These data must 
either be derived adding more error to the calculation or substituted with commanded 
rates. It is important that all known biases be removed from sensor data and calculated 
input errors since the Kalman filter estimator is based on the assumption of uncorrelated 
zero-mean Gaussian noise. Even if all input parameters were known exactly the discrete 


modeling of the spacecraft dynamics introduces some error. 


1. External Control and Disturbance Torques 


In the dynamic gyro, known externally applied moments are integrated in the 
system angular momentum calculations. These mclude control moments other than those 
imparted by momentum exchange devices as well as modeled disturbance torques. Since 
control torques are normally of significant magnitudes, it is essential that they be 
modeled correctly. Moments from magnetic torque rods depend on the imperfectly 
controlled magnetic dipole and the earth’s magnetic field strength. The magnetic dipole 
is provided by torquer current measurement. The earth’s magnetic field must either be 
modeled or measured with magnetometer but is not precisely known. Reaction jet 
moments are almost impossible to model accurately. This often result in degraded 


pointing and tracking during firings operations. 


The effects of external disturbance moments depends on spacecraft configuration 


and orbital profile. Those disturbances that have significant effects on vehicle dynamics 
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should be modeled whenever possible to increase the accuracy of the attitude control 
system. These types of errors are non-zero mean in the short term and therefore can only 
be corrected for with sensor measurement updates. An extended option for The 
Aerospace Corporation’s Pseudo Gyro includes a torque bias estimator to reduce the 


effects of unknown external disturbances. 


For the Bifocal Relay Mirror satellite, the gravity gradient torque is the most 
significant disturbance and can be modeled as an input to increase the accuracy of the 
dynamic rate calculation. The model is well understood and is generated from the inertia 
matrix which is already required by the dynamic gyro and vehicle orientation with 
respect to the gravity vector which can be determined from the estimated attitude and 
ephemeris data. The gravity gradient model used by the dynamic gyro is equivalent to 


the subsystem shown in Figure 13. 


2. Reaction Wheel Relative Momentum 


Instead of integrating the torques produced by momentum exchange devices their 
relative momentum effects are used directly in the dynamic rate calculation. The relative 
momentum of each reaction wheel is given by its orientation within the spacecraft, the 
component inertia of its spinning disk and the wheel spin rate. The imperfect sensor 
measurements from the reaction wheel tachometers introduce errors in system 
momentum calculation. Relative orientation angles of reaction wheels are fixed and 
errors can be corrected through calibration. Orientations of control moment gyros, 
however, are variable. Since these devices usually carry much more momentum small 
gimbal resolver errors can have a significant impact on total system momentum 
calculations. In this simulation, time varying artificial alignment errors are applied to the 
reaction wheel momentum measurements to observe these effects. The error corrupted 
wheel momentum measurement subsystem implemented in the SIMULINK model is 


shown in Figure 32. 
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Figure 32. SIMULINK Subsystem Diagram: Error Corrupted Reaction Wheel 
Momentum Measurement 


3. Moment of Inertia Calculations 


The calculation of the total spacecraft inertia matrix is accomplished by Equation 
(3.12) using the SIMULINK subsystem illustrated in Figure 7. Fixed component 
moments of inertia and masses of the primary and secondary bodies are assumed to be 
known as well as the relative positions of their centers of mass from the system mass 
center. Imperfect knowledge of these parameters introduces errors in the angular rate 
calculation. Errors generated from rotating spacecraft components are not constant in the 


spacecraft body frame. 


The inertia matrix also depends on internal sensor input from position encoders or 
potentiometers for relative angular orientation of appendages. The model of the 
potentiometer that measures the relative angle of the receive telescope includes 
quantization effects and additive noise. If appendage relative motion is slow or 
component moments of inertia are small, it may not be necessary to update the system 
inertia dyadic at the bandwidth of the attitude processor. A trigger is added to the 
SIMULINK subsystem that calculates the inertia matrix (Figure 7), so that the affects of 


update rate can be evaluated. 


4. Appendage Relative Momentum 


Knowledge of appendage relative momentums has direct bearing on the system 
momentum and therefore the dynamic rate calculation. The relative momentum of the 


Bifocal Relay Mirror satellite’s secondary body is calculated from Equation (3.19). It 
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depends on the knowledge of the fixed secondary body inertia matrix and the relative 
angular rate. The attitude simulation model assumes that there is no directed 
measurement of the relative angular rate. The rate must therefore be derived from 
potentiometer measurements of the relative orientation about the axis of rotation. This 
sensor may have a minimum discernable incremental angle and noise corruption. Also, 
since the rate is derived from position measurements it exhibits increased noise and time 
lag. The simulation diagram that produces the error corrupted relative angle and rate for 


the spacecraft secondary body is shown in Figure 33. 
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Figure 33. SIMULINK Subsystem Diagram: Appendage Relative Angle and Rate 
Measurement 


Appendages of significant inertia or relative rates like the receive telescope of the 
Bifocal Relay Mirror satellite should be controlled with smooth gimbal drive motors in 
order to minimize nonlinear relative pointing errors. Ideally, all moveable appendages 


would have an associated rate sensor for each axis of rotation. 


The other option for approximating appendage relative angular rates is to use the 
commanded drive input. In some systems, appendage controllers can provide a smooth 
relative rate through the drive motor, which significantly enhances angular momentum 
tracking. Drive actuators with considerably erratic friction effects will cause errors in the 
attitude control system during slew operations. Short term transients may introduce 


significant settling times for error corrections. 


The subsystem of the attitude determination simulation that controls all of the 
internal sensor measurements and input parameter calculations for the dynamic gyro is 


shown in Figure 34. 
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Figure 34. SIMULINK Subsystem Diagram: Dynamic Gyro Inputs 


D. CALCULATED ANGULAR RATE ERROR CORRECTION 


The error state extended Kalman filter is designed to track rate bias errors that are 
relatively constant in the spacecraft body frame. This is a good approximation for 
properly functioning gyro based systems. For the filter to remain effective when used 
with dynamic gyro, the bias in the output rate as seen in body coordinates due to the 
various error sources must be small and exhibit a bandwidth below the measurement 
update rate. Since it is the case that the common error sources are not zero-mean, their 
effects on the rate output must be relatively constant so that the momentum correction 
supplied by the Kalman filter applies over the update interval. Transient error spiking 
can require multiple star tracker measurement updates to correct. If the internal sensors 


provide inconsistent measurements, dynamic gyro angular rate will be degraded. 
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VII. RESULTS 


Simulation results demonstrate that multi-body spacecraft attitude control without 
the use of rate gyroscopes can be performed using the attitude and angular rate estimation 
scheme proposed in this thesis. The performance of dynamic gyro based attitude 
determination system is compared with a similar gyro based system using the Bifocal 
Relay Mirror attitude simulation. Simulation input parameters are varied to analyze the 
effects of major error sources on the dynamic gyro model. Also evaluated are the effects 
of star tracker accuracy and measurement update rates on the attitude determination 


system. 


A. BASELINE SIMULATION 


A full set of MATLAB plots are presented in this section for baseline analysis of 
the dynamic gyro based attitude determination and control scheme developed in this 
thesis. This set of results validates the potential effectiveness of the proposed attitude 
and angular rate estimating scheme used for multi-body spacecraft control. It also 
provides a common reference for analysis and comparison with subsequent simulation 
results. For other simulations, only selected plots that are required for analysis of results 


will be shown. 


1. Simulation Input Parameters 


Table 3 shows the inputs parameters that are held constant for each simulation run 
used to obtain results. These input parameters are set in the MATLAB script file that 
calls the SIMULINK attitude simulation. The MATLAB code file is included as 

















Appendix B. 
SIMULATION PARAMETERS (SIMULINK) 
Simulation Time Period 500 sec 
Attitude Determination Bandwidth 20 Hz 
SIMULINK Solver Method ode5 (Dormand-Prince) 
Solver Fixed Step Size 0.05 sec 


ORBIT PARAMETERS (Circular Orbit) 





715 km 


65 


COMMANDED MANEUVER PROFILE 
Inertial Attitude Quaternions See Figure 35 


Body Axes Angular Rates See Figure 36 
Secondary Body Relative Angle See Figure 37 
DISTURBANCE MOMENTS 

Gravity Gradient Modeled 

Secular (Magnitude) le-4 Nm 

Periodic (Magnitude, Period) 4e-4 Nm, [400,500,600] sec 
Disturbance Effect on System Angular See Figure 38. 

Momentum 























MOMENTS OF INERTIA MATRICIES 
Primary Body [2997 .28025,- 
3.9331,118.2824 
-3.9331,3164.18285,1.1230 
118.2824, 1.1230,881.82105] 





ken? 
Secondary Body [1721.07340,-0.0116,- 
7.8530 
-0.0116,1559.85414,- 
12.5463 
-7.8530,- 
12.5463,182.89962] kgm” 


CENTER OF MASS OFFSET FROM SPACECRAFT C.M. 


Primary Body [0.558354158, 
3.91788e-4, 
0.15226902] m 


Secondary Body [-1.302113918, 
-9.13673e-4, 
0.355100092] m 

WHEEL CONTROL LAWS 

Quaternion Error Gains (Kq) [3000,7000,4500] 

Angular Rate Error Gains (Kw) [1000,2000, 1000] 

Control Law Delay for Initial 30 sec 

Determination Errors 


GYROSCOPE CHARACTERISTICS 
































Static Rate Biases le-4*[-1,1.5,1] rad/sec 
Rate Noise Variance le-8 
Acceleration Noise Variance (Rate le-12 
Random Walk) 
Quaternion Errors (q1,q2,q3) 
Angular Rate Errors [-0.001,0.001,0.002] rad/sec 
KALMAN FILTER INITIALIZATION 
State Estimate | [0,0,0,0,0,0] 
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State Error Covariance Matrix 5e2*[100,0,0,0,0,0 
0,100,0,0,0,0 
0,0,100,0,0,0 
0,0,0,1,0,0,0 
0,0,0,0,1,0,0 
0,0,0,0,0,1,0 
0,0,0,0,0,0, 1] 

Tracker | (x-rotation, y-rotation) 135 deg, 30 deg 

Tracker 2 (x-rotation, y-rotation) 135 deg, -30 deg 

Tracker 3 (x-rotation, y-rotation) 180 deg, 0 deg 

REACTION WHEEL PARAMETERS (4 WHEEL PYRAMID 
CONSTELLATION) 


Number of Wheels Operating 
Constellation Angle to xy-Plane 45 deg 




















Constellation Torque Saturation Limit 1 Nm 
Constellation Torque Rate limits 10 Nm/sec 
Table 3. Simulation Input Parameters 


2. Command Attitude Profile 


The 500 second maneuvering profile chosen for the Bifocal Relay Mirror attitude 
simulation analysis is illustrated in Figures 35, 36 and 37. This profile resembles the 
maneuver required to maintain transmit and receive telescope pointing control during an 
overhead operational pass to conduct laser relay operations. The majority of the 
maneuver is performed in the spacecraft pitch axis, q), as both telescopes orient to point 
at fixed ground sites. Less significant motion is required in the spacecraft roll and yaw 
axes in order to ensure that the relative axis of rotation of the receive telescope is 
correctly oriented during the tracking maneuver. Based on ground site separation 
distance and orbital altitude, the largest relative angle require between the telescopes is 
about 30 degrees during a near overhead pass between the uplink and downlink ground 


sites. 
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Commanded Quaternions 









































Figure 35. Baseline Commanded Attitude Profile 



































Figure 36. Baseline Commanded Angular Rate Profile 
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Figure 37. Baseline Commanded Relative Angle to Receive Telescope Profile 


The magnitude of the total spacecraft angular momentum during the maneuvering 
profile is plotted in Figure 38. The spacecraft attitude is completely controlled by 
momentum exchange with the reaction wheels. Therefore, the changes to the angular 
momentum profile are attributable the external disturbance moments. The gravity 
gradient disturbance is modeled and the other disturbances are fixed so the angular 


momentum profile remains essentially common in all simulation runs. 
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Figure 38. Total Spacecraft Angular Momentum Profile 
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3. Attitude Determination System Performance Results 


In this baseline simulation, updates from randomly selected star trackers are 
provided to the dynamic gyro and the attitude propagator via the Kalman filter at two 
second intervals. The star tracker H and V measurements are corrupted with noise 
variance of 1x10*. In order to observe the performance of the attitude determination 
system without updates, a 200 second star gap is simulated starting 100 seconds into the 
run. As a worst-case analysis, this star gap occurs during the peak maneuvering time of 
the satellite including the rotation of the secondary body. Attempt is made to tailor the 
plant and measurement error covariance matricies used in the Kalman filter. The attitude 


determination system is initiated with the errors given in Table 3. 


The accuracy of the angular rate calculation is entirely dependent upon the ability 
of the dynamic gyro to track the total spacecraft angular momentum. The error in the 
magnitude of the total system momentum compared to the simulated actual momentum is 
shown in Figure 39. The steady state momentum error is held within 0.07 Nms with 
consistent star tracker data but builds to 0.35 Nms after 200 seconds without stars. No 


external disturbance torques are modeled as dynamic gyro inputs in this simulation run. 
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Figure 39. Baseline Dynamic Gyro Angular Momentum Error 


After star tracker measurements are processed, the Kalman filter provides a 


momentum correction to the dynamic gyro determined from the spacecraft inertia and the 
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rate bias error state estimate. The attitude quaternion is also updated and the states are 


reset to zero. Figure 40 shows the magnitude of the error states determined by the 


Kalman filter. No updates are provided during the 200 second star gap. 
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Figure 40. 


x10 State Attitude Errors 


















0 50 100 150 200 250 300 350 400 450 500 
Time (sec) 
x10 State Rate "Bias" Errors 








0 50 100 150 200 250 300 350 400 450 500 
Time (sec) 





Baseline Kalman Filter Attitude and Rate Bias Errors 


Quaternion errors in the attitude determination system are plotted in Figure 41. 


Steady state attitude errors are maintained within 3x10~ during periods of continuous star 


coverage but are increased an order of magnitude by the end of the 200 second star gap. 


After the star gap, the attitude error build up is quickly removed through measurement 


corrections. 
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Figure 41. Baseline Estimated Attitude Quaternion Error 











Figure 42 illustrates the nature of the angular rate estimation error produced by 
the dynamic gyro based determination system. The error generated in the roll axis, w,, is 
the most significant because it is aligned with the axis of rotation of the secondary body. 
The relative rate is derived from imperfect potentiometer measurements. The rate errors 
increase during receive telescope motion due to the potentiometer quatization effect 
which produces the broken pattern of noise. The yaw axis is most coupled dynamically 
to roll axis and exhibits similar errors at less magnitude. The build up of the bias in the 
rate error is hard to perceive among the noise but the effects are evident in the quaternion 


error plot. 
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Figure 42. Baseline Estimated Angular Rate Error 


4. Attitude Control System Performance Results 


Attitude control is accomplished by three reaction wheels selected from a pyramid 
constellation of four. Nonlinear response of these control actuators is simulated by 
saturation and rate limiting the output torques as shown in Figure 11. Wheel torque 
commands are generated by combining feed forward and error based control laws. The 
control laws implementation is delayed at initiation to allow the attitude determination 
system to converge. High gain control laws dominate the wheel torque response. Figures 
43 and 44 show the torque and momentum response of the three operating reaction 


wheels. 
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Figure 43. Reaction Wheel Control Torques 
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Figure 44. Reaction Wheel Angular Momentum 


Attitude control performance is dominated by attitude determination errors. The 
control system is designed to limit steady state errors in the attitude quaternions to less 
than 2x10’ with ideal attitude and angular rate knowledge. With the control laws 
applied, the attitude control response resembles the errors of the attitude determination 
system in steady state. Figure 45 shows the attitude quaternion error response for the 
baseline simulation run. Damped corrections back to steady state after the star gap are 


much slower than those realized by the attitude determination system. 
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Figure 45. Baseline Control Attitude Quaternion Error 


The error in the controlled spacecraft angular rate is shown in Figure 46. Since 
most of the attitude determination noise is removed by the controller, the uncorrected rate 
error bias due to the star gap is evident. Although the errors remain small during the star 


gap, the bias drives the attitude quaternion error build up. 
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Figure 46. Baseline Control Angular Rate Error 





B. DYNAMIC GYRO VS GYRO PERFORMANCE 


Results presented in this section compare a gyro based attitude determination 


system with the dynamic gyro based system. Three orthogonally mounted gyroscopes 


D 


are modeled to provide simulated spacecraft body angular rate measurements to the 
attitude propagator when the gyro option is selected. The gyro model includes static bias, 
rate noise, and rate random walk supplied by integration of random noise. Gyroscope 
characteristics used for the simulation results are listed in Table 3. A bias error tracking 
system is also modeled in parallel with the gyro generated rates to accept bias corrections 
provided by the Kalman filter. Direct comparison of the dynamic gyro and gyro based 
attitude determination systems is accomplished using identical input parameters to obtain 
simulation results. The SIMULINK subsystem for gyro rate measurement and bias error 


tracking is illustrated in Figure 47. 


NON-IDEAL GYRO MEASUREMENT 
AND BIAS ERROR TRACKING 
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Figure 47. SIMULINK Subsystem Diagram: Gyro Measurement and Bias Error 
Tracker. 
1. Continuous Star Tracker Coverage 


With precise attitude updates from star trackers there is only little noticeable 
difference between the dynamic gyro and gyro based attitude determination systems. 
Figure 48 provides a side-by-side comparison of simulation results obtained with a 


continuous star update interval of 2 seconds. 
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Results: Dynamic Gyro vs. Gyro with Continuous Star Coverage 


Figure 48. 


Although the attitude quaternion errors show little difference between the two 
systems the rate errors from the gyros are noticeably more accurate and better resemble 
white noise. Since the characteristics of the rate errors are different between the two 
systems the plant error covariance matricies used in the Kalman filter are chosen 
differently. Table 4 shows the plant error covariance matricies used by the dynamic gyro 
and gyro based systems. These values were determined through tuning over several 


simulation runs. 
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Plant Error Covariance Matrix (Q) 
Dynamic Gyro 
0 





Table 4. Plant Error Covariance Matricies 
2. Gapped Star Tracker Coverage 


Figure 49 shows the comparison of the dynamic gyro and gyro based attitude 
determination systems when a 200 second gap in star tracker coverage is encountered 
during maneuvering operations. Very little difference is evident in the rate error plots 
because the bias build up during the star gap is so small compared to the noise in the 
error. The effect of the bias, however, shows in the attitude quaternion error plots. At 
the end of the star gap the error in the dynamic gyro based system is about five times that 


of the gyro based system. 
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Figure 49. Results: Dynamic Gyro vs. Gyro with Gapped Star Coverage 


c. DYNAMIC GYRO PLANT ERROR ANALYSIS 


In this section the effects of major error sources on the performance of the 
dynamic gyro are analyzed. Figure 50 shows the attitude determination plots for the 
baseline simulation run. Dynamic gyro input parameters are varied and results are 


compared with these plots. 
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Figure 50. Baseline Attitude Determination Performance Results 


1. Disturbance Torque Modeling 


The accuracy of the dynamic gyro depends on knowledge of external disturbance 
moments and internal spacecraft momentum. For the baseline simulation run no external 
disturbances were modeled in the dynamic gyro. This leads to the relatively large drift in 
attitude quaternions when uncorrected during the star gaps. If the spacecraft’s attitude 
with respect to the orbital reference frame can be determined the gravity gradient 
disturbance torque can be modeled. For the Bifocal Relay Mirror satellite this 
disturbance has the greatest effect. Figure 51 illustrates the increase in performance of 
the dynamic gyro based attitude determination system when the gravity gradient moment 
is modeled. At the end of the star gap the attitude errors are comparable to the gyro 
based system. The total system angular momentum error is much smaller during the star 


gap and improvement can even be seen during times of continuous star coverage. 
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Figure 51. Gravity Gradient Disturbance Modeled in the Dynamic Gyro 
2. Rotation Axis Alignment Error 


Alignment errors of momentum exchange control devices and slewing 
appendages have direct effects on the dynamic gyro momentum error. Alignment errors 
of fixed reaction wheels are easily compensated for by calibration but control moment 
gyros have directionally variant momentum vectors with respect to the spacecraft body. 
The results plots shown in Figure 52 are produced when a periodic alignment error on the 
net momentum of the reaction wheels with a magnitude of approximately 0.5 degrees is 
added to the simulation. A significant increase in attitude error is developed during the 
star gap. The dynamic gyro does not track the system angular momentum as well even 


during continuous star coverage. 
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Figure 52. Reaction Wheel Alignment Error Effects on the Dynamic Gyro 
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3. Potentiometer Quantization 


The patterns shown in angular rate error plots for the dynamic gyro do not appear 
as white noise during appendage motion because of quantization effects in the model of 
the potentiometer that measures the relative orientation. The simulated quantization level 
is 0.01 degrees in the baseline simulation. Since the effect alternates direction of error 
the attitude quaternion errors are not affected significantly. Figure 53 shows simulation 
performance results when the quantization level is decreased an order of magnitude. The 


angular rate error looks more like white noise and the bias build up during the star gap 


can be seen. 
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Figure 53. Reduced Potentiometer Quantization Effect on Dynamic Gyro 
Performance 


4. Moment of Inertia Update Frequency 


The inertia matrix for spacecraft with relatively small or slowly moving 
appendages does not change quickly. In these cases processing power may be saved by 
performing the system inertia calculation at a lower bandwidth than the dynamic gyro. In 
the Bifocal Relay Mirror satellite the secondary body is large so even small slew rates 
cause significant change in the system inertia matrix. Figure 54 shows the effects of 
decreasing the inertia calculation frequency from the 20 Hz rate of the dynamic gyro to 
once every 10 seconds. The quaternion error profile is significantly altered but the 
magnitude of the error is only slightly increased. The momentum error in the dynamic 
gyro takes longer to correct after the star gap since the Kalman filter attempts to correct 


for all errors as if they were due to spacecraft momentum. 
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Figure 54. 10 Hz Moment of Inertia Calculation 


D. STAR TRACKER MEASUREMENT ERROR ANALYSIS 


The dynamic gyro based attitude determination system is highly reliant on 
Kalman filter updates based on star tracker measurements. The quality of the Kalman 
filter corrections depend on the accuracy of the trackers, the discrete measurement 


interval and number of trackers providing measurements. 


if Star Tracker Accuracy 


The results shown in Figure 55 are from simulations using random selection at 2 
second intervals. The plot on the left is generated using a 1 sigma variance of 1x10“ in 
the horizontal and vertical measurements of the trackers. To generate the plot on the 
right the variance is reduced by four times. The direct effect on attitude determination 


performance is apparent. 
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Figure 55. 


2. Update Interval 


Figure 56 shows a comparison of attitude determination performance using 2 and 
6 second star measurement intervals. Shorter update intervals result in quicker 


convergence and less steady state error. 
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Figure 56. Star Tracker Update Interval: 2 Seconds [left] vs 6 Seconds [right] 








3. Star Tracker Selection 


Because the star tracker measurements provide useful information along only two 
axes at least two trackers must be active to maintain pointing control. Figure 57 shows 
the comparison of an attitude determination simulation with updates spread evenly 
between three trackers and one that uses two trackers with 95% of the updates coming 
from the same sensor. There is no apparent degradation when the updates are spread 


evenly between the two operational trackers. 
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Figure 57. Star Tracker Selection Comparison: Even Distribution [left] vs One 
Tracker Favored [right] 


85 


THIS PAGE INTENTIONALLY LEFT BLANK 


86 


VIII. CONCLUSIONS 


Dynamic modeling provides an imperfect but operative means of estimating 
multi-body spacecraft angular rates when mechanical gyro data are not available. An 
attitude and angular rate estimation scheme is developed in this thesis that integrates the 
dynamic gyro concept with an error state extended Kalman filter estimator that provides 
precise attitude updates from star tracker measurements. Results indicate that the 


determination system provides effective estimates for performing attitude control. 


A. SUMMARY 


The attitude determination system is incorporated into a multi-body spacecraft 
attitude simulation for evaluation and analysis. Simulation is an extremely valuable 
analysis tool for understanding the effects of error sources on system performance. It is 
also a suitable mechanism for comparison with gyro based determination systems since a 


mechanical gyro model can easily be inserted. 


The effects of the primary error sources on the dynamic gyro performance are 
investigated through simulation. It is shown that the corrections provided by a star 
tracker based Kalman filter make the system robust to measurement and parameter 
knowledge error sources. Significant improvement in attitude determination performance 
is realized when disturbance torques are modeled. The other primary error sources 
include the alignment of momentum exchange control devices and relative angle and rate 
knowledge in large or quickly slewing appendages. Error effects are amplified during 


star gaps when no corrections to the dynamic model are available. 


The software implemented dynamic gyro essentially emulates the functions of a 
set hardware gyroscopes. In a spacecraft where the mechanical gyros have failed or 
become too erratic to be corrected by the Kalman filter the dynamic gyro may be a viable 
replacement. Operated in tandem with mechanical gyros, either system can provide 


redundant inertial angular velocity for improved attitude control. 
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This attitude determination concept is ideally suited to a spacecraft designed 
specifically for its implementation with precise internal sensors and mechanisms to 
monitor spacecraft parameters and integrated external torque estimation modeling. In 
these systems the dynamic gyro can increase the lifetime and reliability of the spacecraft 


while reducing the power requirements. 


B. RECOMMENDATIONS 


A strong recommendation for future work is simulation model tailoring to real 
spacecraft hardware. The multi-body dynamics can be expanded to model an actual 
satellite’s mass and inertia characteristics. Additionally, sensor parameters can be 
modeled to match existing sensor error specs. With these modifications the simulation 
model can be used conduct analysis and predict performance when the dynamic gyro 
software is implemented on board the spacecraft. Future work may include using the 
simulation to compare the predicted dynamic gyro based attitude determination 


performance to actual gyro based performance recorded with telemetry data. 


There are several improvements that can be made to the attitude determination 
simulation model. It could easily be modified to provide the capability to mix gyroscope 
and dynamic gyro data for redundant information. A parity matrix developed from the 
pseudo-inverse concept can be generated to account for system observability in the over 
determined case. This would allow performance analysis with selected gyro failures and 
further aid in the evaluation of its utility as a back-up attitude determination scheme. 
Another simulation improvement would be the incorporation of higher order dynamic 
effects into the model including center of mass offsets and flexibility modes for the 


appendage couplings. 


Finally, improvements to the attitude determination scheme can be developed 
including a torque error estimator as suggested by The Aerospace Corporation [Ref. 8]. 
The better the dynamics are modeled in the processing software the more effective the 


dynamic gyro becomes as a replacement for the hardware gyroscope. 
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APPENDIX A: KALMAN FILTER BACKGROUND 


“The Kalman filter combines all available measurement data, plus prior 
knowledge about the system and measuring devices, to produce an estimate of the desired 
variables in such a manner that the error is minimized statistically” [Ref. 2]. An 
understanding of the Kalman filter requires some background in the theory of probability 


of random variables and processes [Ref. 1]. 


A. PROBABILITY 


The probability of an event, e, represents a possible outcome of a random 
experiment and is written Pr(e). A random variable, X, can be thought of as a function of 
the outcomes of some random experiment. The manner of specifying the probability with 
which different values are taken by a random variable is the probability distribution 


function, F(x) 


F(x)=Pr(X < x) (A.1) 
This is a function over the range of possible values that shows the probability with which 
the random variable takes on a value at or less than the value of the range. Its derivative 
is the probability density function, f(x), 
dF(x) 
dx 


This function identifies the likelihood of a random variable assuming a particular value in 


(A.2) 





f(x)= 


its range of possible values. Over the range of all possible values a characteristic of any 


probability distribution or density function is 


Flos) | f(u)du =1 (A.3) 


A joint probability density function can be defined for multiple variables. For two 


random variables the joint density is given by 


0°F, (x,y) 


a (A.4) 


f (x,y) oa 
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B. EXPECTATION 


The expectation of a random variable is defined as the sum of all values the 
random variable may take, each weighted by the probability with which the value is 
taken. 


EIXI= xf(x)dx (A.5) 


-0co 


which is also called the mean value or first moment of X. This is a precisely defined 
number toward which the average of a number of observations of the random variable X 
tends. Since a function of a random variable is itself a random variable, the expectation 


of a function of X, E[g(X)], can be expressed as 


Elg(X)I=| 2(x)f(x)dx (A.6) 


An important statistical parameter descriptive of X is its mean squared value defined by 


E[X?] = | xf(x)dx (A.7) 


—oo 


which is also called the second moment of X. The root-mean-squared (rms) value of X is 
J E[X’] . The variance of a random variable is defined as the mean squared deviation of 


the variable from its mean denoted by 6”. 
6° = | (K-E[X]) f(x)dx = E[X* -E[XP (A.8) 


For zero mean random variables the variance is simply E[X*]. The square root of the 
variance, 6, is called the standard deviation of the random variable. 

A very important concept is that of statistical correlation between random 
variables. A partial indication of the degree to which one variable is related to another is 


given by the cross covariance, which is the expectation of the product of the deviations of 


the two variables from their means, 
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E[(X-E[X])(Y-E[Y])]= } | (X-ELX])(y-ELY Df, (%y)dydx=E[XY]-E[X]E[Y] 


348 
348 


(A.9) 


For a vector of random variables a symmetric covariance matrix can be defined 
where the diagonal elements are the individual variances of the vector components and 
the off diagonals are given by the cross covariances between the corresponding vector 
components. The cross covariance, normalized by the standard deviations of X and Y, is 
called the correlation coefficient. 


_ E[XY]-E[X]JE[Y] 


— (A.10) 
OxOy 


This is a measure of the degree of linear dependence between variables X and Y. If they 


are completely independent fi is zero. 


C. LEAST-SQUARES ESTIMATION CONCEPT 


The optimality condition used in the Kalman filter is the minimization of 
weighted least-squares error. The least-squares minimization problem involves a set of 
measurements, y, which are linearly related to the vector of state variables, x , by the 
expression 

y=Hx+v (A.11) 
where V is an unknown vector of actual measurement errors with zero mean. The error 
we seek to minimize is the measurement residual, Z, given by 

Z=y-Hx (A.12) 


where X is the estimate of the actual state vector. Since the sum of the squares of a 


vector are generated by the inner product, the cost function, J, is 


J=(y-H&) (y-HX) (A.13) 


Minimization is obtained by differentiating and setting the result to zero 
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Ca, (A.14) 


Ox 


and ensuring that the Hessian is positive semidefinite 


a°J 


ae 20 (A.15) 











The resulting least-squares solution for the state estimate is 


=(H"H) H'y (A.16) 
A weighted least-squares solution can be used when not all components of the 


measurement residual are treated equally 


%=(H™WH) | H"wy (A.17) 
where W is a positive semidefinite matrix relating scale factors between components. 
Expanding this solution to multiple estimates over time assumes all measurements are 
used together in a batch processing scheme. Every time a new measurement becomes 
available it is appended to the measurement vector and the estimated state vector includes 
estimates corresponding to each of the accumulated measurements. The Kalman filter is 
based on recursive processing where each measurement is used sequentially to generate 
an optimal estimate of the current state without recomputing estimates of all previous 
states. For this technique all previous information is embodied in the prior estimate and 


state covariance matrix. 


D. STATE ESTIMATE AND COVARIANCE 


The Kalman filter estimation algorithm maintains the first two statistical moments 


of the estimated state. The estimated state, X , IS a vector of random variables whose 
mean (first moment) is the actual state vector, x. The error in the state estimate is 


defined as 


x Ux —-X (A.18) 














and is thus assumed to be zero mean. ‘The covariance of the state estimate error, 


designated P, is given by 
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P=E[XX"] (A.19) 
This covariance matrix associated with the state estimate error (second moment) provides 


a Statistical measure of the uncertainty in x and the correlation between the errors of its 


components. For a state vector of two variables the state estimate error is given by 


i= * (A.20) 
X 
and has the covariance matrix 
P=E Be a s Elsa. BE i (A.21) 
XK 5, Ky E[x,x,]  E[x;] 


The diagonal elements of the state covariance matrix are the mean square errors in 
the knowledge of the state variables while the off diagonal elements are indicators of 


cross correlation between elements of the estimated error. 


E. STOCHASTIC LINEAR DYNAMIC MODEL 


The Kalman filter requires representation of system dynamics in a linearized 
state-space form, a linear measurement model, and assumed characteristics of process 
and measurement noise. For a continuous linear system the general state-space model 


and measurement equations are given by 


X(t)=F(t)X+G(t)w(t) 
y(Q=H(Hx(H)+v(t) 


where w(t) and v(t) are random vectors representing the unmodeled disturbance inputs 


(A.22) 


and measurement errors. These vectors are treated as unbiased (zero mean) white noise. 
Note that if one of these vectors is known to have a nonzero bias the mean can be 
augmented onto the state vector creating a new state space model with unbiased random 


error. The covariance matricies associated with the process disturbance, w(t), and the 


measurement noise, V(t), are 
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Q=Elv(tyvt)"] 


. (A.23) 
R=E[w(t)w(t) J 
The equivalent discrete representation is 
XO, E+ A x 
QPHATAW x (A.24) 


Yi =H, X FV, 
Here the subscript k represents values at time { and subscript k+1 represents values at the 
next discrete time step t,,. Note that for time invariant systems the discrete state 


transition matrix, 6, , is related to the continuous formulation by 


Oraen (A.25) 
which depends only on the systems dynamics matrix, F, and the discrete time interval. If 
the discrete interval is kept short enough the relation holds as an approximation. The 
state transition matrix allows the calculation of state vector at the next discrete time step 


in the absence of forcing functions. It obeys the differential equation 


= = F(t), (A.26) 





F. PROPAGATION OF ERRORS 


The recursive Kalman filter requires the propagation estimate and_ error 
covariance based on system dynamics. In the discrete implementation the error in the 


current estimate given by 

X, =%,-X, (A.27) 
has the covariance matrix representing the uncertainty in the estimate 

PSE x] (A.28) 
The best estimate of the future state, Ka , 1s given by 

Xpu= Ox (A.29) 
The error in the new estimate 

Xp = OX, — A Wy (A.30) 


has the expected value 
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E[x,,,]= 6, Elx, + A E[w, ]=0 (A.31) 
since X, and W, are assumed to be unbiased. Thus the predicted estimate remains 
unbiased. Note that if a deterministic input is added to the dynamic system model then 
the identical quantity is added to both the actual and estimated state leaving the 


estimation error unchanged. It can be shown that the associated state error covariance 


matrix is given by 


Pu = FLX Xi] 5, P.O; a A Q rN (A.32) 
where Q is the covariance of the random system disturbance. It is evident from this 
equation that the size of disturbance directly impacts the uncertainty in the state estimate. 
The system dynamic stability reflected in the state transition matrix also effects the 
uncertainty. The covariance of neutrally stable or unstable systems will grow unbounded 
over time in the absence of measurement corrections. The propagation of the state 
estimate and its error covariance matrix is the prediction step of the Kalman filter 


algorithm. 


G. MEASUREMENT UPDATES 


The correction step of the Kalman filter algorithm incorporates measurement 
updates to refine the state estimate and error covariance. This step is executed only when 
a measurement becomes available. When a measurement is taken we use the superscripts 
‘-° and ‘+’ to denote values at a particular time before and after incorporation of the 
measurement correction. Given a prior estimate of the system state at time { we seek to 


update our estimate based on the measurement y, in a linear, recursive form 


x, =K,X, + Ky, (A.33) 
where K, and K, are matricies, as yet unspecified, that determine the relative weighting 


of the prior estimate and current measurement. The error in this estimate can be shown to 


be 
xt =(K, +K HA, -1)%,+K,x, +K,¥, (A.34) 


Since xX, and V, are unbiased 
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E[x; ]=0 


if (A.35) 

E[V,]=0 
this estimate can only remain unbiased for all given states if 

K, =I-K,H, (A.36) 
With this requirement the estimator takes the form 

Xf =x, +K, (¥, Hix; ) (A.37) 
and has corresponding error 

xf =(I-K,H,)x,+K,¥, (A.38) 


The state error covariance must also be updated. 


P= Ex x, | 


= E{(d—-K,H,)x,[x, (@-K,H,)" + ¥,K, ]+K,¥,[x," d-K,H,)’ +9, K; ]} 


(A.39) 

Using the definitions 

P, = E[x,x,] (A.40) 
and 

R, =E[V,%, ] (A.41) 
and noting that the measurement errors are uncorrelated 

E[x, ¥, ]=E[¥,x,']=0 (A.42) 
The covariance can be simplified to 

P, = (1-K,H,)P, (I- K,H,)' + K,R,Ky (A.43) 


The criterion for choosing the optimal K, is to minimize a weighted scalar sum of 


the diagonal elements of the error covariance matrix, P,”. Thus the cost function is 


J, = EES Xi] (A.44) 


where S is any positive semidefinite matrix. Choosing S = I, 
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J, =trace[P, ] (A.45) 
which is equivalent to minimizing the length of the estimation error vector. To find the 
value of K which yields a minimum, it is necessary to take the derivative of J} with 


respect to K, and set it equal to zero. The result is 


—2(I— K, H, )P-H; + 2K,R, =0 (A.46) 
Solving for Kx, 
K, =P, H{(H,P, HT" (A.47) 


which is referred to as the Kalman gain matrix. Since the Hessian of J, is positive 
semidefinite so K does indeed produce a minimum. Substituting into the equation for 


the updated error covariance gives 


Pt =P, — P-H?[H,P,H! +R, J"H,P; = (I-K,H, )P; (A.48) 
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APPENDIX B: MATLAB CODE 


A. INPUT PARAMETERS AND SIMULATION CALL 


%Bifocal Relay Mirror Spacecraft Attitude Control System 
% Input parameters and control options for attitude simulation model 


clear 
tic 


stoptime=500; %Simulation stoptime 
dt=.05; %fixed step size and computer time step 


m1=2267.6059; %Mass of X-mit telescope (kg) 
m2=972.3628; Z%oMass of RCV telescope (kg) 


%Orbital parameters (circular) 

Re=6378.1363; %Earth radius (km) 

mu=3.986004415e5; % gravitational constant (km‘43/s‘*2) 

h=715; %orbit altitude (km) 

r=Re+th; %orbit radius (km) 

we=7.2921158533e-5; %Earth rotation rate (r/s) 

w0=sqrt(mu/r*3); %orbital rate based on circular orbit/altitude (1/s) 

nu0=0; %initial true anomaly (r) 

u0=0; initial angle of the magnetic dipole normal to the 
%line of ascending nodes (r) 

el=11*pi/180; %magnetic dipole tilt angle 

inc=40*pi/180; Y%inclination 

Km=7.943e15; %magnetic constant (Nm‘2/a) 

Bm=Km/(r*1000)43; Ymagnetic field strength (N/am) 


%Disturbance torques 

Mds=1e-4*[1;1;1]; %secular disturbance torques (Nm) 
Mdp=4e-4*[1,1,1]; %periodic disturbances torque magnitudes (Nm) 
pdper=[400,500,600]; %periodic disturbance torque periods (s) 
pdfreq=2*pi./pdper; Yoperiodic disturbance torque frequencies (r/s) 
pdphase=[pi/4,pi/2,pi]; Yoperiodic disturbance torque phase (r) 


%Magnetic dumping control 

magoo=0; %magnetic dumping on/off control (O=off, 1=on) 
magsat=180; %maximum torque rod output 

magk=5e5; Yomagnetic torquer gain 


% Appendage (RCV scope) motion relative to body x-axis (sinusoidal) 
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bper=300; %period of oscillation (s) 

bstart=.25*bper; Yostart motion (s) 

bmotion=1*bper; %length of motion (s) 

bstop=bstart+bmotion; %stop motion (s) 

bfreq=2*pi/bper; Yofrequency of motion (rad/s) 

bamp=15*pi/180; %amplitude/half maximum angular offset (rad) 

bi=0; %initial angular offset (rad) 

bdi=0; %initial relative angular rate (rad/s) 

bacc=bamp*bfreq*2; % maximum realtive angular acceleration (rad/s’2) 
bnv=le-5“2; %relative angular noise variance from potentiometer (rad) 
bns=0; %initial seed of potentiometer noise 

bq=.01*pi/180; Yoquantization of potentiometer readings (rad) 


%Command profile inertial to body accelerations (sinusoidal) 
wdcfreq=2*pi*[1/800; 1/500; 1/600]; %frequency of acceleration profile variation 
wdc=[.05*wdcfreq(1)*2;.5*wdcfreq(2)*2;-.1*wdcfreq(3)%2]; 

Yamplitude of acceleration profile 
wdcphase=[pi/8;-pi/12;0]; %phase of acceleration profile 
cpoo=1; %turn commanded profile on/off (O=off, 1=on) 


11=[2997.28025,-3.9331,118.2824; 

-3.9331,3164.18285,1.1230; 

118.2824,1.1230,881.82105]; 

Ymoment of inertia matrix of primary body (XMIT scope) about its cm 
12=[1721.07340,-.0116,-7.8530; 

-.0116,1559.85414,-12.5463; 

-7.8530,-12.5463,182.89962]; 

Ymoment of inertia matrix of appendage (RCV scope) about its cm 


dx 1=.558354158; %Distance of X-mit telescope cm from S/C cm in x dir (m) 
dy 1=3.91788e-4; %Distance of X-mit telescope cm from S/C cm in y dir (m) 
dz1=-.15226902; %Distance of X-mit telescope cm from S/C cm in z dir (m) 
dx2=-1.302113918; %Distance of RCV telescope cm from S/C cm in x dir (m) 
dy2=-9.13673e-4; %Distance of RCV telescope cm from S/C cm in y dir (m) 
dz2=.355100092; %Distance of RCV telescope cm from S/C cm in z dir (m) 


IC1=I+m1*[dy14%2+dz1%2,dx1*dy1,dx1*dz1; 
dx1*dy1,dx142+dz142,dy1*dz1; 
dx1*dz1,dy1*dz1,dx142+dy1%2]; 
Ymoment of inertia matrix of primary body (XMIT scope) about S/C cm 
roti=[1,0,0;0,cos(bi),sin(bi);0,-sin(bi),cos(bi) |; 
%rotation matrix corresponding to initial relative angular position 
%of appendage (RCV scope) 
%x-axis rotation of magnitude bi 
Ic2i=roti'*I2*roti; Zinitial rotated moment of inertia matrix of appendage 
%(RCV scope) 
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I2m=m2*[dy242+dz242,dx2*dy2,dx2*dz2; 
dx2*dy2,dx2%2+dz2%2,dy2*dz2; 
dx2*dz2,dy2*dz2,dx2%2+dy2%2]; 
Ymoment of inertia of appendage (RCV scope) about S/C C.M. 
%due to mass offset 
IC2i=roti'*I2*roti+I2m; %initial moment of inertia matrix of appendage 
%(RCV scope) about S/C cm 
li=IC14+IC2i; %initial S/C moment of inertia about S/C cm 


% Wheel control law gains 
You=-Kqqe-Kwwe=-Kq(qce*q)-Kw(w-wc) 

cldel=30; %delay control laws for attitude determination 
Kq1=3000; 

Kq2=7000; 

Kq3=4500; 

Kw1=1000; 

Kw2=2000; 

Kw3=1000; 


%Initial body with respect to inertial quaternions 
qli=0;q2i=0;q3i=0; 
qi=[q1i,q2i,q3i,sqrt(1-qli*2-q2i*2-q3i*2)]; 


%Initial body with respect to orbit quaternions 
qloi=0;q20i=0;q30i=0; 
qoi=[q10i,q201,q301,sqrt(1-q1loi*2-q20i1*2-q30i%2)]; 


%Initial body with respect to inertial direction cosine matrix 
demi(1,1)=qi(1)*2-qi(2)*2-qi(3)42+qi(4)%2; 
demi(1,2)=2*(qi(1)*qi(2)+qi3)*qi(4)); 
demi(1,3)=2*(qi(1)*qi(3)-qi(2)*qi(4)); 
demi(2,1)=2*(qi(1)*qi(2)-qi(3)*qi(4)); 
demi(2,2)=-qi(1)*2+qi(2)*2-qi(3)42+qi(4)%2; 
demi(2,3)=2*(qi(2)*qi(3)+qi(1)*qi(4)); 
demi(3,1)=2*(qi(1)*qi(3)+qi(2)*qi(4)); 
demi(3,2)=2*(qi(2)*qi(3)-qi(1)*qi(4)); 
demi(3,3)=-qi(1)*2-qi(2)*2+qi(3)%2+qi(4)%2; 
Ci=[dcemi(1,1);dcemi(1,2);dcemi(1,3); 

demi(2,1);demi(2,2);dcemi(2,3); 

demi(3,1);dcemi(3,2);dcemi(3,3)]; 

% Direction cosine component vector 
Ai=[demi(1,1),dcmi(1,2),dcmi(1,3); 

demi(2,1),dcmi(2,2),dcmi(2,3); 

demi(3,1),dcmi(3,2),dcmi(3,3)]; 

% Direction cosine/attitude matrix 
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%Initialize S/C momentum with wi at orbital rate 
wi=[0;-w0;0]; %Initial angular rate of S/C primary body due to orbital rate 
Hi=li*wi; %initial angular momentum of S/C due to orbital rate 


%Gyro characteristics 
grsb=le-4*[-1,1.5,1]; %arbitrary gyro static bias 
grv=le-9*[1,1,1]; variance of gyro rate noise 
grm=[0,1,2]; %initial seed of gyro rate noise 
gasb=[0,0,0]; %gyro acceleration static bias (zeros) 
gav=le-12*[1,1,1]; %variance of gyro acceleration noise 
gan=[3,4,5]; %initial seed of gyro acceleration noise 
grrwi=[0,0,0]; Yoinitial rate random walk (zeros) 
wbi=[0;0;0]; %initial gyro bias correction 
gx=0*pi/180; %body to gyro(IRU) x-axis rotation 
gy=0*pi/180; %body to gyro(IRU) y-axis rotation 
G=[cos(gy),sin(gy)*sin(gx),-sin(gy)*cos(gx); 

0,cos(gx),sin(gx); 

sin(gy),-cos(gy)*sin(gx),cos(gy)*cos(gx)]; 

%Gyro alignment matrix - body to IRU 


% Dynamic Gyro input options 

dg=1; %choose Dynamic Gyro or regular gyros (1=DG, 0=gyros) 

dgge=0; %calculate gravity gradient torque for DG (1=yes(requires or=1), 0=no) 
dgsd=0; %secular disturbance torque known for DG (1=yes, 0=no) 

dgpd=0; %periodic disturbance torque known for DG (1=yes, 0=no) 

moiu=1; Y%steps per periodic moment of inertia calculation (1=every step) 

or=0; %perform orbital reference calculations 


% Attitude determination initialization with error 
ge li=q1i+.008;qe2i=q2i+.012;qe3i=q3i-.008; 
qei=[ge li,qe21,qe31,sqrt(1-ge1i*2-qe2i*2-qe3i*2)]; 
dcmei(1,1)=qei(1)*2-qei(2)*2-qei(3)*2+qei(4)’2; 
dcmei(1,2)=2*(qei(1)*qei(2)+qei(3)*qei(4)); 
dcmei(1,3)=2*(qei(1)*qei(3)-qei(2)*qei(4)); 
dcmei(2, 1)=2*(qei(1)*qei(2)-qei(3)*qei(4)); 
dcmei(2,2)=-qei(1)42+qei(2)*2-qgei(3)*2+qei(4)*2; 
dcmei(2,3)=2*(qei(2)*qei(3)+qei(1)*qei(4)); 
dcmei(3, 1 )=2*(qei(1)*qei(3)+qei(2)*qei(4)); 
dcmei(3,2)=2*(qei(2)*qei(3)-qei(1)*qei(4)); 
dcmei(3,3)=qei(1)*2-qei(2)*2+qei(3)2+qei(4)*2; 
Cei=[dcmei(1,1);dcemei(1,2); 
dcemei(1,3);dcmei(2, 1);dcmei(2,2);dcmei(2,3); 
dcemei(3,1);dcmei(3,2);dcmei(3,3)]; 
Aei=[dcmei(1,1),dcmei(1,2),dcmei(1,3); 
dcemei(2,1),dcmei(2,2),dcmei(2,3); 
dcmei(3,1),dcmei(3,2),dcmei(3,3)]; 
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geloi=qloi+.02;qe20i=q20i+.03;qe30i=q30i-.02; 

qeoi=[qel o1,qe201,qe30i,sqrt(1-geloi*2-ge20i*2-ge301%2)]; 
wei=witle-3*[-151;2]; 

Hei=hi* wei; 


% Initialize Kalman filter 
xi=0*[ones(3,1)*eps* 1e10;ones(3,1)*eps*1e8]; %initial state vector 
Pi=[100*eye(3),zeros(3);zeros(3),eye(3)]*500; Yinitial covariance matrix 
if dg== 
Qn=.5; 
else %plant noise covariance constant DG vs gyro 
Qn=.03; 
end 
Q=[100*eye(3),zeros(3);zeros(3),eye(3)]*Qn; Yplant noise covariance 
%small Q = good plant, big Q = bad plant 
R=10000*eye(2); %measurement noise covariance (10000) 
E=[1,0,0;0,1,0]; %choose horizontal and vertical components only 


%Star tracker parameters 
corr=2* 1/dt; Ycomputer steps between stars 
gapstart=100; %start gap in star tracker data (s) 
stargap=200; %length of star gap (s) 
gapstop=gapstart+stargap; Y%end gap in star tracker data (s) 
tlx=135*pi/180; %body to star tracker one x-axis rotation 
tly=30*pi/180; %body to star tracker one y-axis rotation 
T1=[cos(tly),sin(tly)*sin(t1x),-sin(tl y)*cos(t1x); 
0,cos(t1x),sin(t1x); 
sin(tly),-cos(tly)*sin(t1x),cos(tly)*cos(t1x)]; 
%body to star tracker one rotation matrix 
t2x=135*pi/180; Ybody to star tracker two x-axis rotation 
t2y=-30*pi/180; %body to star tracker two y-axis rotation 
T2=[cos(t2y),sin(t2y)*sin(t2x),-sin(t2y)*cos(t2x); 
0,cos(t2x),sin(t2x); 
sin(t2y),-cos(t2y)*sin(t2x),cos(t2y)*cos(t2x)]; 
%body to star tracker two rotation matrix 
t3x=180*pi/180; %body to star tracker two x-axis rotation 
t3y=0*pi/180; %body to star tracker two y-axis rotation 
T3=[cos(t3y),sin(t3y)*sin(t3x),-sin(t3y)*cos(t3x); 
0,cos(t3x),sin(t3x); 
sin(t3y),-cos(t3y)*sin(t3x),cos(t3y)*cos(t3x)]; 
%body to star tracker two rotation matrix 
mixT3=1; %choose whether or not to mix in star tracker 3 data (l=yes, 0=no) 
T3th=.3; Y%choose threshold for random mix of star tracker 3 data 
%(0=50%, higher=less, lower=more) 
T12b=0; %choose to bias random selection between tracker | and 2 
%(0=50/50, lower=favor1, higher=favor2) 
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stnb=[0,0]; Z unknown star tracker bias error 
stnv=le-4%2*[1,1]; %variance of star tracker noise 
stns=[0,1]; %initial seed of star tracker noise 


%3 of 4 reaction wheels in pyramid constellation 

rwsat=1*[1,1,1]; %constellation wheel torque saturation limits (Nm) 

rwrl=10*[1,1,1]; Zconstellation wheel torque rate limiter (Nm/s) 

hwi=-12*[bdi;0;0]; %initialize wheels to cancel appendage momentum 

rwa=45*pi/180; %angle of reaction wheels to x-y plane 

RW=inv([-cos(rwa)*cos(pi/4),-cos(twa)*cos(pi/4),cos(rwa)*cos(pi/4); 
-cos(rwa)*cos(pi/4),cos(rwa)*cos(pi/4),cos(rwa)*cos(pi/4); 
-sin(rwa),-sin(rwa),-sin(rwa)]); Z%body to wheel transform matrix 

hwsi=RW*hwi; %distribute initial wheel momentum 

rwnm=[0;0;0]; %mean of wheel noise 

rwnv=le-242*[1;1;1]; %variance of wheel noise 

rwns=[7;8;9]; %initial seed of wheel noise 

rwev=le-8*[1,1,1]; %variance of wheel alignment walk error noise 

rwes=[10,11,12]; %initial seed of gyro acceleration noise 

rwam=1*5e-3*[1,2,-1]; %reaction wheel constellation alignment error magnitudes 

rwaf=2*pi*[100,50,20]; %reaction wheel constellation alignment error frequencies 

rwap=[pi/4,pi/3,pi/2]; Yoreaction wheel constellation alignment error phase 


sim(‘acs_sim509') %Call to SIMULINK simulation 


%Call ploting files for simulation results analysis 
profileplots %command profile 

ADplots509 %attitude determination performance 
ACplots509 %attitude control performance 


toc 


B. SIMULATION RESULTS PLOTS 


1. Commanded Profile Plots 


%Laser Realy Spacecraft command profile plots 


figure(1) 
clf 
plot(tout,qc) 
title(Commanded Quaternions') 
legend(‘q1','q2','q3','q4') 
xlabel('Time (sec)') 
grid on 
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figure(2) 

clf 

plot(tout,wc) 

title((Commanded Angular Rates') 
legend(‘wl','w2','w3') 
xlabel('Time (sec)') 
ylabel(‘rad/s') 

grid on 

figure(3) 

clf 

plot(tout,b* 180/pi) 

xlabel('Time (sec)') 

ylabel('‘Deg') 

title('Relative Angle to Recieve Telescope (alpha)’) 
grid on 

figure(4) 

clf 

Hmesqrt(H(:, 1).42+H(:,2).42+H(:,3).%2); 
plot(tout,Hm) 

title(‘System Angular Momentum’) 
xlabel('Time (sec)') 

ylabel(‘Nms') 

grid on 


2 Attitude Determination Performance Results 


%Bifocal Realy Mirror satellite attitude determination plots 
%analyze performance of dynamic gyro and Kalman filter 


%DG momentum error or gyro bias 

if dg== 

figure(5) 

clf 

Hme=esqrt(H(:, 1).42+H(:,2).42+H(:,3).%2); 
Hmdg=saqrt(Hdg(:,1).42+Hdg(:,2).42+Hdg(:,3).%2); 
plot(tout,Hm-Hmdg) 

title’ Angular Momentum Error’) 
xlabel('Time (sec)') 

ylabel(‘Nms') 

axis([0,stoptime,-.4,.4]) 

grid on 

else 

figure(5) 

clf 

plot(tout,wb) 

title(‘'Gyro Biases’) 
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legend(‘wb1','wb2','wb3’) 
xlabel('Time (sec)') 
ylabel(‘rad/s') 

grid on 

end 


%KF states 

figure(6) 

clf 

subplot(211),plot(tout,x(:, 1:3)) 
title(‘State Attitude Errors’) 
legend(‘xal','xa2','xa3') 
xlabel('Time (sec)') 
ylabel(‘rad') 
axis([0,stoptime,-.0008,.0008]) 
grid on 
subplot(212),plot(tout,x(:,4:6)) 
title('State Rate "Bias" Errors’) 
legend(‘xb1',’xb2','xb3') 
xlabel('Time (sec)') 
ylabel(‘rad/s') 
axis([0,stoptime,-.00008,.00008]) 
grid on 


%AD errors 

figure(7) 

clf 

plot(tout,(q-qkf)) 
title(Quaternion Errors’) 
legend(‘q1','q2','q3','q4') 
xlabel('Time (sec)') 
axis([0,stoptime,-3e-3 ,3e-3]) 
grid on 

figure(8) 

clf 

plot(tout,w-wkf) 

title('Rate Errors’) 
legend(‘wl','w2','w3') 
xlabel('Time (sec)') 
ylabel(‘rad/s') 
axis([0,stoptime,-1.5e-3,1.5e-3]) 
grid on 


3. Attitude Control Performance Results 


% Bifocal Realy Mirror satellite attitude control plots 
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Yanalyze commanded versus controled attitude 


figure(9) 

clf 

plot(tout, Mws) 

title("‘Wheel Torques') 
xlabel('Time (sec)') 
ylabel('Nm') 

legend(‘Wheel 1',"Wheel 2','Wheel 3’) 
grid on 

figure(10) 

clf 

plot(tout,hws) 

title(Wheel Angular Momentums') 
xlabel('Time (sec)') 
ylabel(‘Nms') 

legend(‘Wheel 1',"Wheel 2','Wheel 3’) 
grid on 

figure(1 1) 

clf 

plot(tout,(q-qc)) 
title(‘Quaternion Error’) 
legend(‘q1','q2','q3','q4') 
xlabel('Time (sec)') 
axis([0,stoptime,-4e-3 ,4e-3]) 
grid on 

figure(12) 

clf 

plot(tout,w-wc) 

title(‘Rate Error’) 
legend(‘wl','w2','w3') 
xlabel('Time (sec)') 
ylabel(‘rad/s') 
axis([0,stoptime,-1.5e-3,1.5e-3]) 
grid on 
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